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Abstract 

This  study  examines  the  performance  and  limitations  of  a  heuristic  cooperative 
control  (CC)  surveillance  algorithm  for  multiple  unmanned  aerial  vehicles  (UAVs)  un¬ 
der  both  simulation  and  demonstration.  The  algorithm  generates  Dubin’s  based  paths 
and  provides  velocity  feedback  to  accomplish  simultaneous  arrival  onto  a  surveillance 
orbit  around  the  target  and  maintains  position  while  orbiting.  The  performance  was 
tested  under  multiple  wind  conditions  in  simulation  and  actual  winds  during  flight 
testing.  Both  position  accuracy  and  target  visibility  were  examined. 

The  analysis  covers  three  major  topics:  development  of  a  closed  loop  model  for 
a  new  airframe  at  AFIT  for  simulation  purposes,  development  of  the  CC  algorithm 
that  interfaces  with  Procerus  Technologies’  Kestrel  Autopilot,  and  achievable  system 
performance  analysis.  The  model  assumes  first-order  responses  to  roll,  pitch,  and 
airspeed  commands  using  time  constants  pulled  from  actual  flight  test  data.  The 
CC  algorithm  has  two  modes:  one  that  generates  commands  to  multiple  UAVs  for 
simultaneous  arrival  to  a  surveillance  orbit,  and  one  that  maintains  equal  angular 
spacing  about  the  orbit.  In  addition  to  positional  performance  metrics,  percentage  of 
target  in- view  time  was  also  measured  based  on  the  UAV’s  side  camera  held  of  view 
(FOV).  Simulation  tested  both  modes  under  wind  conditions  of  0%,  10%,  25%,  and 
50%  of  the  nominal  airspeed  ( Vnom ). 

Results  showed  that  the  algorithm  maintained  UAV  position  with  winds  25% 
of  Vnom,  but  instabilities  appeared  at  50%  where  large  overshoots  appeared  on  the 
downwind  side  of  the  orbit.  Target  visibility  was  most  impacted  by  crosstrack  errors 
that  steadily  grew  with  increasing  winds.  Roll  of  the  UAV  showed  the  greatest  impact 
on  the  FOV  due  to  its  coupling  effect  with  crosstrack  error.  Overall  target  in-view 
time  also  improved  with  increasing  numbers  of  UAVs  for  all  wind  conditions. 
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Multiple  Unmanned  Aerial  Vehicles 


I.  Introduction 

1.1  Background 

The  Air  Force  is  utilizing  Unmanned  Aerial  Vehicles  (UAVs)  at  an  ever  in¬ 
creasing  pace.  Small  autonomous  vehicles  have  sparked  great  interest  in  the  military 
by  providing  an  inexpensive  system  that  increases  capabilities  and  prevents  placing 
personnel  in  dangerous  situations.  Autonomous  platforms  that  fly  have  a  unique  ap¬ 
peal.  They  can  traverse  large  distances  quickly  and  provide  a  “bird’s  eye  view”  of 
the  battlespace.  Utilizing  multiple  vehicles  enhances  mission  accomplishment  with 
redundancy,  robustness,  and  increased  coverage  when  compared  to  a  single  platform. 
This  research  explores  operating  multiple  UAVs  for  surveillance. 

Multi-UAV  surveillance  holds  many  advantages  over  the  other  surveillance  op¬ 
tions  in  terms  of  proximity  (close  or  far),  speed  (fast  or  slow),  responsiveness,  cost, 
and  overall  personnel  risk.  Manned  surveillance  is  close  and  responsive  but  slow  and 
places  personnel  at  risk,  traditional  aircraft  are  fast  and  reasonably  close  but  are  ex¬ 
pensive  and  also  place  personnel  at  risk.  Space  surveillance  allows  access  to  denied 
areas  but  is  very  expensive,  limited  by  the  orbit  for  timing  and  placement,  and  far 
from  the  target.  Multi-UAV  surveillance  can  reduce  or  remove  personnel  risk,  be  close 
to  the  target,  provide  persistence  over  the  target,  and  can  cost  very  little  compared 
to  manned  aircraft  and  space  options.  Advances  in  the  miniaturization  of  electron¬ 
ics  aided  this  interest  in  UAVs.  As  cost  and  size  decreased,  capability  increased  for 
surveillance  and  autonomous  technology.  Consequently,  research  and  development 
blossomed  in  both  academia  and  the  aerospace  industry. 

The  Air  Force  Institute  of  Technology  (AFIT)  has  conducted  many  UAV  re¬ 
search  projects,  and  vigorously  continues  to  this  day.  The  research  conducted  herein 
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continues  an  ongoing  project  at  AFIT  that  focuses  on  utilizing  UAVs  for  surveillance 
and  target  engagement  missions.  To  fill  the  void  in  data  for  small  aircraft,  AFIT’s  re¬ 
search  began  in  2006  with  the  work  of  Nidal  Jodeh  [8]  modelling  a  9.16  foot  wingspan 
radio  controlled  model  airplane,  the  Sig  Rascal,  retrofitted  with  an  autopilot.  The 
stability  and  payload  capacity  of  this  airframe  made  it  ideal  for  UAV  research.  Since 
the  Sig  Rascal  became  the  primary  demonstration  aircraft,  Jodeh’s  model  became 
the  base  for  many  following  projects.  One  UAV  application  focused  on  tracking  and 
engaging  a  moving  target  with  on-board  video.  At  AFIT,  this  application  became 
known  as  the  “Fleeting  Target  Program.”  The  problem  was  broken  up  into  creating  a 
path  to  the  target  in  real  time  (Pathmaker)  [19],  flying  the  vehicle  to  the  target  using 
video  feedback  (Cursor  On  Target)  [20],  and  integrating  the  hardware  and  software 
into  a  usable  package  (Fleeting  Target  Controller)  [17].  The  research  described  herein 
is  the  next  iteration  of  the  Fleeting  Target  Program. 

For  the  current  effort,  the  emphasis  of  research  shifted  away  from  target  detec¬ 
tion  and  engagement  when  Air  Force  Research  Lab  (AFRL)  received  an  urgent  need 
request  from  the  warfighter  to  develop  a  route  surveillance  capability.  AFRL  was 
tasked  to  deliver  a  prototype  system  capable  of  monitoring  many  miles  of  road  and 
revisiting  any  point  at  fixed  intervals  [2],  The  proposed  system  consists  of  multiple 
UAVs  with  day  and  night  sensors,  a  ground  station  with  semi-autonomous  control  of 
the  UAVs,  and  an  anomaly  detection  system.  The  primary  purpose  was  to  surveil 
routes  ahead  of  convoys  to  minimize  risk  to  transportation  operations.  This  need 
became  the  primary  drive  for  this  research. 

To  aid  this  research,  AFRL  provided  AFIT  with  six  new  airframes,  the  Bat¬ 
tlefield  Air  Targeting  Camera  Autonomous  Micro  UAVs  (or  BATCAMs  for  short). 
This  airframe  differed  greatly  from  Jodeh’s  airframe:  The  BATCAM  wingspan  was 
21”  vs  Jodeh’s  9.1  ft,  the  BATCAM’s  propulsion  was  electric  vs.  Jodeh’s  gas  engine, 
cameras  were  body  fixed  vs.  gimballed,  and  control  surfaces  were  a  V-tail  with  no 
ailerons  vs.  the  traditional  aileron/rudder /elevator  configuration.  The  old  Sig  Rascal 
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model  no  longer  applied.  By  changing  the  demonstration  aircraft,  research  using  the 
BATCAM  dictated  development  of  a  new  model. 

The  scope  of  this  work  develops  and  assesses  two  main  items:  a  model  for  the 
BATCAM,  and  an  algorithm  that  controls  multiple  BATCAMs  for  surveillance. 

1.2  Problem  Statement 

This  research  focuses  on  quantifying  the  abilities  and  limitations  of  a  group  of 
UAVs  to  monitor  a  fixed  target.  The  primary  focus  is  not  target  detection,  but  the 
control  algorithm  for  the  group  that  accomplishes  persistent  visual  contact  of  the 
target.  Quantifying  the  abilities  and  limitations  of  the  UAV  group  are  done  from  two 
perspectives:  simulation,  and  flight  demonstration. 

The  simulation  portion  requires  construction  of  a  model  for  the  BATCAM, 
then  controlling  this  model  with  an  algorithm  to  accomplish  surveillance.  Using 
flight  test  data,  this  research  constructs  a  representative  closed-loop  model  for  the 
BATCAM/autopilot  system.  The  intent  of  this  model  is  to  capture  the  major  handling 
characteristics  of  the  system,  then  utilize  an  algorithm  that  focuses  on  UAV  placement 
with  respect  to  one  another  in-flight  to  surveil  a  target. 

The  flight  demonstration  portion  utilizes  only  the  algorithm  to  command  the 
UAV,  replacing  the  model  with  an  actual  BATCAM  and  autopilot.  The  algorithm 
will  analyze  the  current  state  of  multiple  UAVs  and  create  a  command  set  to  place 
all  UAVs  for  surveillance.  The  algorithm  accomplishes  two  distinct  tasks:  to  place  all 
UAVs  into  the  surveillance  orbit,  and  to  maintain  the  surveillance  orbit. 

The  goal  of  this  research  is  to  command  multiple  UAVs  in  real  time  to  approach 
and  maintain  an  orbit  about  a  fixed  target.  The  effects  of  wind  on  performance 
as  well  as  the  ability  of  the  UAV  group  to  reconfigure  when  individual  UAVs  are 
added/deleted  will  be  determined. 
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1.3  Assumptions  and  Limitations 

To  keep  a  concise  focus,  certain  bounds  and  simplifications  must  be  made.  Re¬ 
search  of  this  nature  is  highly  complex,  and  exhaustive  treatment  of  a  subject  can 
quickly  exceed  the  time  and  resources  available.  The  limitations/bounds  will  apply 
to  both  simulation  and  demonstration,  but  the  assumptions  may  only  apply  to  the 
mo  del/ simulation.  If  one  of  the  items  applies  to  simulation  or  demonstration  only, 
the  list  will  specify  the  applicability. 

The  research  bounds  are: 

•  The  “changing”  conditions  will  be  limited  to  wind  scenarios,  and  adding  and 
removing  a  UAV  under  wind. 

•  The  maximum  number  of  UAVs  will  be  four. 

•  All  UAVs  will  be  the  same. 

The  following  lists  the  overall  assumptions.  Chapter  III  will  elaborate  on  the  reasoning 
behind  these  items. 

•  Each  UAV  closed-loop  model  will  behave  like  a  first  order  system  (for  simulation 
portion  only). 

•  Wind  vectors  will  be  only  in  the  horizontal  plane  (simulation  only). 

•  The  ground  will  be  assumed  planar  and  flat  for  sensor  footprint  projections. 

•  Communication  is  available  to  all  UAVs  at  all  times  (simulation  only). 

•  The  flat  earth  model  will  be  assumed  “inertial”  ignoring  the  rotation  and  cur¬ 
vature  of  the  earth. 

The  intent  of  these  assumptions  is  to  make  simplifications  that  are  reasonably  accurate 
to  actuality  but  have  the  net  benefit  of  decreasing  complexity.  The  real  time  aspect  of 
the  algorithm  relies  on  accurate  yet  timely  information  flow  to  and  from  UAVs.  The 
objective  of  this  research  is  to  create  a  heuristic  real-time  algorithm  that  is  robust 
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enough  to  handle  environmental  effects.  Comparisons  to  truly  optimal  solutions  will 
be  accomplished  in  future  research. 

1.4  Preview 

Chapter  II  presents  past  research  from  two  areas:  background  research  that 
this  thesis  builds  upon,  and  related  research  that  illustrates  the  different  ways  to 
control  and  optimize  multiple  UAVs  under  different  scenarios.  Chapter  III  presents 
the  equipment  used  for  flight  demonstration,  the  modelling  approach  for  simulation, 
and  the  heuristic  algorithm.  Chapter  IV  presents  the  performance  results  of  the 
heuristic  algorithm  under  simulation  and  flight  demonstration.  Chapter  V  concludes 
this  research  and  makes  recommendations  for  future  research. 
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II.  Literature  Review 


2. 1  Overview 

This  survey  draws  from  two  areas:  background  research,  and  related  research. 
The  background  research  pulls  together  material  that  aids  this  thesis.  Whereas  the 
related  research  explores  how  various  researchers  approached  the  multi-UAV  coop¬ 
erative  control  problem.  With  a  very  different  airframe  from  past  Fleeting  Target 
research,  the  first  portion  of  background  research  aids  in  developing  a  new  model. 
Jodeh  [8]  creates  a  framework  for  modelling  small  vehicles,  then  Kingston  [9]  from 
AFRL  uses  a  simple  closed-loop  model  for  the  route  surveillance.  The  next  area  of 
background  research  is  path  planning.  From  the  knowledge  base  here  at  AFIT,  the 
works  of  Zollars  and  Terning  illustrate  two  differing  approaches:  one  optimal  and  one 
heuristic.  Rysdyk  [16]  from  the  University  of  Washington  also  developed  a  nice  closed 
form  solution  for  keeping  a  target  in  the  field  of  view  (FOV)  while  orbiting.  The  re¬ 
lated  research  section  surveys  some  applications  of  multiple  UAVs  used  cooperatively. 
Secondly,  Beard  [4]  uses  the  broad  area  search  technique  to  compare  optimal  verse 
sub-optimal  solutions.  Next,  MIT’s  Richards  and  Bellingham  [15]  take  the  scenario  of 
completing  X  tasks  using  UAVs,  and  increase  the  complexity  to  include  differing  UAV 
capabilities,  time  constraints  between  steps,  and  no  fly  zones.  Ending  this  section  is  a 
brief  overview  of  a  powerful  simulation  tool  called  MultiUAV2.  This  product  creates 
a  realistic  environment  to  test  algorithms  and  highlights  not  only  optimal  path  plan¬ 
ning,  but  also  key  considerations  for  any  fielded  system  like  probability  of  detection 
and  inter-UAV  communication. 

2.2  Background  Research 

Creating  a  good  mathematical  description  of  the  BATCAM  airframe  is  not  a 
minor  feat.  Both  the  approach  and  the  aerodynamic  properties  must  be  thorough  for 
generating  a  representative  model.  The  combined  works  of  Jodeh  [8]  and  Kingston  [9] 
provide  a  solid  base  to  build  from. 
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2.2.1  Aircraft  Modelling-  Jodeh.  Jodeh’s  [8]  research  created  one  of  the  first 
models  of  a  micro  air  vehicle  (MAV)  at  AFIT.  He  found  that  extensive  research  existed 
for  modelling  large  airframes,  but  high  fidelity  small  vehicle  models  seemed  rare.  Not 
surprisingly,  he  also  found  no  research  on  UAV  handling/stability  characteristics  or 
standards  for  handling/stability,  both  important  to  accomplish  quality  surveillance 
missions.  The  need  for  a  modelling  approach  for  small  vehicles  seemed  apparent.  The 
Sig  Rascal  airframe  by  Tower  Hobbies  was  modelled.  This  is  a  stable  airframe  with 
the  ability  to  handle  small  to  medium  payloads.  A  Piccolo  autopilot  gave  the  airframe 
autonomy. 

This  modelling  approach  did  not  use  the  wind  tunnel,  but  used  the  United 
States  Air  Force  (USAF)  Stability  and  Control  Digital  Datcom  software.  This  com¬ 
puter  program  was  written  for  the  Air  Force  under  contract  by  McDonnell  Douglas 
Astronautics  Company.  This  program  has  the  ability  to  output  Lift,  Drag,  Moment, 
and  Stability  Derivatives  (and  many  other  items)  when  given  desired  flight  condi¬ 
tions,  attitudes,  and  physical  geometry.  With  this  output,  Jodeh  created  a  6  Degree 
of  Freedom  (DOF)  non-linear  model  in  Matlab/Simulink.  The  model  was  compared 
to  both  flight  testing  and  Hardware-in-the-loop  (H1TL)  simulations. 

To  create  a  good  comparison  between  the  three  sets  of  data,  he  used  an  elevator 
command  which  induced  both  a  phugoid  and  short  period  motion.  Comparing  the 
flight  test  data  to  the  6  DOF  model,  the  actual  (flight  test)  oscillation  had  a  period 
of  13s  where  the  model  predicted  10s.  In  the  phugoid  excitation,  the  model  matched 
in  period,  but  the  amplitude  differed:  flight  test  values  exceeded  the  model  in  both 
pitch  rate  and  airspeed,  but  model  amplitude  values  exceeded  test  values  in  altitude 
oscillations.  A  comparison  of  Hardware  in  the  loop  (HITL)  to  test  showed  larger 
frequencies  and  less  damping  for  HITL.  All  difference  were  within  a  factor  of  two, 
but  model  values  and  HITL  values  differed  from  test  values  up  to  25%  in  period  and 
50%  in  amplitude.  For  these  types  of  measurements,  being  within  a  factor  of  two  is 
reasonable  and  indicates  that  the  model  is  fairly  close.  Small  errors  in  the  parameters 
of  a  non-linear  model  can  manifest  like  differences  shown  above.  Since  MAVs  have 
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small  masses  and  moments  in  inertia,  errors  in  aerodynamic  forces  and  moments  make 
state  errors  more  pronounced. 

Jodeh’s  research  lays  a  good  foundation  to  follow  for  modelling  small  vehicles 
for  this  research.  Not  only  does  it  provide  a  fairly  complete  list  of  items  to  model, 
but  also  illustrates  that  a  model  for  a  small  vehicle  can  differ  from  reality  within  an 
order  of  magnitude  from  expected  values. 

2.2.2  UAV  Closed  Loop  Model/Route  Surveillance-Kingston.  Members  of 
AFRL  Vehicles  Directorate  [9]  presented  both  a  simple  UAV  model  and  an  algo¬ 
rithm  that  is  inherently  decentralized,  convergent  to  optimal  behavior  in  finite  time, 
accounts  for  communication  range  limitations,  and  allows  for  changing  perimeters. 

Each  UAV  uses  this  simplified  model  derived  from  the  6  DOF  model.  The  model 
assumes  constant  altitude  and  constant  airspeed.  The  autopilot  has  also  been  tuned 
so  that  the  closed-loop  vehicle  behaves  like  a  first  order  system.  The  equations  of 
motion  become 

PN  =  VaCOS'lp  +  WN 

Pe  =  vasimp  +  wE 

ip  =  Pjr tamp  (2.1) 

*  a 

Va  =  kv(VaC  -  Va) 

p  =  k^((pc  -  0) 

where  and  p^  are  position;  -0,  0,  V  are  yaw  angle,  roll  angle,  and  airspeed;  [wjv, 
We]  is  the  wind  vector;  k ’s  are  the  first  order  parameters,  and  c  denotes  commanded 
quantities.  This  is  a  nice  simple  approach  to  capture  the  behavior  of  an  airframe/au¬ 
topilot  closed-loop  system.  This  equation  will  be  modified  to  create  the  BATCAM 
model  in  Chapter  III. 

Kingston  goes  on  to  use  this  model  to  surveil  a  stretch  of  road.  For  N  UAVs, 
the  road  is  divided  into  N  equal  segments.  With  all  vehicles  on  a  single  line,  they 
travel  back  and  forth  along  their  assigned  segment  and  meet  at  certain  times.  It  is 


possible  for  2  (or  more)  vehicles  to  occupy  the  same  position  on  the  line.  To  handle 
communications,  the  analysis  requires  that  UAVs  physically  meet  along  the  line  to 
exchange  information.  The  base  location  is  located  at  one  end  of  the  line  (x=0), 
so  when  a  UAV  reaches  the  end,  information  reaches  the  base.  As  intuition  would 
hold,  by  evenly  spacing  all  vehicles,  all  information  can  reach  any  part  of  the  line  in 
minimum  time. 

Each  UAV  patrols  the  ith  segment  back  and  forth.  When  the  timing  works 
perfectly,  the  two  UAVs  will  briefly  meet  at  the  endpoints  and  exchange  information. 
If  the  first  UAV  meets  the  second  UAV  early  (or  within  its  segment),  the  second  UAV 
will  turn  around  and  be  escorted  back  to  the  end  of  the  segment,  then  turn  around. 
If  the  second  UAV  is  late,  the  first  UAV  will  proceed  into  the  other  UAV’s  segment 
until  they  meet.  At  that  time,  the  first  UAV  will  turn  around  and  be  escorted  back 
to  the  division  point.  As  they  contact,  they  exchange  information,  both  surveillance 
info  and  4  variables.  These  coordination  variables  are:  Pr  and  P l  are  the  perimeter 
lengths  right  and  left  of  the  UAV;  Nr  and  are  the  number  of  UAVs  right  and  left 
of  the  UAV.  This  information  tells  the  UAV  where  it  lies  on  the  total  perimeter  P  (= 
Pi?  +  Pl)  and  what  its  segment  is. 

Let’s  take  the  worst  case  where  the  perimeter  changes  and  the  UAVs  do  not 
know  how  many  total  UAVs  are  patrolling  the  perimeter  P,  and  initial  UAV  positions 
are  random.  If  time  T  is  P/V,  then  according  to  this  research,  all  UAVs  will  reach 
the  “low  latency”  configuration  along  P  in  less  than  5T.  It  will  take  at  most  3T  for  all 
UAVs  to  exchange  information  and  get  the  correct  perimeter  and  number  of  vehicles, 
then  2T  more  to  reach  the  configuration  where  all  UAVs  are  patrolling  equal  segments. 

The  algorithm  was  validated  through  flight  testing  using  two  UAVs  with  Kestrel 
Autopilots.  The  largest  discrepancy  between  simulation  and  real  data  occurred  in  the 
turn  around.  The  required  U  turn  at  the  end  of  the  segments  took  longer  than 
predicted.  The  shared-border  position  of  the  two  UAVs  was  approximately  60%  of  P, 
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compared  to  the  predicted  50%.  The  wind  that  day  was  35%  of  the  airspeed  of  the 
UAVs. 

Kingston’s  research  utilizes  a  simple  closed-loop  model  for  a  road  surveillance 
algorithm.  Comparing  this  model  to  flight  test  data,  the  simulation  results  compared 
fairly  closely  to  flight  test  telemetry.  This  model  is  the  base  approach  taken  to  create 
a  model  in  this  research.  The  control  algorithm  is  one  approach  to  decentralize 
the  perimeter  surveillance  problem  with  enough  robustness  to  handle  both  changing 
perimeter  and  number  of  UAVs.  The  major  difference  between  Kingston’s  surveillance 
method  to  this  one  is  the  communication  constraint.  He  accounts  for  the  fact  that  not 
all  UAVs  are  reachable  at  any  given  time,  whereas  this  research  assumes  connectivity 
at  all  times  to  all  vehicles.  To  extend  the  multi-UAV  surveillance  over  large  distances, 
this  aspect  must  be  addressed. 

2.2.3  Path  Planning- Zollars.  Michael  Zollars  [21]  analyzed  optimal  path 
planning  for  a  single  vehicle  to  a  static  point.  The  analysis  assumed  the  aircraft  is  a 
point  mass,  aircraft  airspeed  is  constant,  and  both  wind  heading  and  wind  magnitude 
are  constant.  Important  to  note  is  that  the  wings  were  assumed  level  during  turns, 
so  the  side  sensor  look  angle  did  not  vary.  The  wind  varied  from  all  directions  and 
magnitudes  were  varied  for  windspeed/airspeed  ratios  up  to  0.7  (wind  is  70%  of  UAV’s 
airspeed).  Initial  headings  were  constrained  such  that  the  UAV  always  took  off  into  the 
wind,  and  the  distance  between  the  beginning  UAV  position  and  target  was  constant 
for  all  runs.  Three  problems  were  presented: 

•  Finding  the  optimal  path  minimizing  flight  time  given  an  initial  heading  and 
position,  and  a  final  position  and  heading. 

•  Finding  the  optimal  path  minimizing  flight  time  given  through  an  urban  canyon 
with  the  initial  and  final  conditions  from  above. 

•  Finding  the  optimal  path  to  both  reach  the  target  and  orbit  the  target  mini¬ 
mizing  flight  time  and  maximizing  time  in  view  while  orbiting. 
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The  first  two  problems  yielded  expected  results.  Flight  times  increased  with  increasing 
windspeed/airspeed  ratios.  Minimum  flight  time  occurred  with  a  tail  wind  (and  vice 
versa).  The  interesting  results  came  when  the  UAV  orbited  the  target.  If  the  sensor 
was  fixed  to  the  airframe  and  the  sensor  footprint  was  a  point,  the  optimization 
scheme  revealed  it  was  physically  impossible  to  track  the  target  100%  of  the  time. 
When  the  footprint  was  modelled  as  a  circle  with  diameter  of  64%  of  the  orbit  radius, 
the  target  was  in  view  for  82%  of  the  orbit  flight  path. 

This  result  highlights  that  even  under  constant  wind  conditions  with  many 
simplifying  assumptions,  orbiting  under  constant  altitude  has  limits  to  keep  the  target 
in  view.  The  windspeed/airspeed  ratio  was  0.2  for  these  results. 

For  multiple  UAVs  in  wind,  certain  viewing  directions  may  not  be  able  to  main¬ 
tain  the  target  in  the  field  of  view.  To  offset  this  another  UAV  must  provide  a  different 
look  angle  to  maintain  persistence.  For  this  research,  the  key  result  is  that  multiple 
UAVs  orbiting  a  target  has  limitations  in  certain  directions  under  windy  conditions. 

2.2.4  Path  Planning-  Terning.  Building  upon  Zollars,  Terning  [19]  created  a 
heuristic  approach  to  generate  a  path  to  engage  a  moving  target  in  constant  wind  for 
real  time  use.  This  iterative  approach  uses  the  straight  line  distance  to  the  target  as 
an  initial  guess,  then  uses  the  calculated  intercept  time  for  the  next  iteration  until  the 
time  for  both  projected  paths  (target  and  UAV)  fall  within  a  specified  tolerance.  This 
iteration  was  subject  to  both  targets  changing  direction,  and  also  latency  in  target 
information  flow.  The  solution  (called  Pathmaker)  is  sub-optimal,  but  is  fast  enough 
for  real  time  use. 

The  concept  is  quite  simple  and  is  illustrated  in  Figure  2.1.  The  first  step  finds 
the  straight  line  distance  to  the  moving  target.  The  next  step  accounts  for  a  maximum 
effort  turn  and  finds  a  new  distance.  Using  the  time  of  travel  for  the  new  path,  the 
projected  target  distance  is  recalculated  (Step  3).  The  UAV  refines  the  needed  path  in 
Step  4.  Realize  this  assumes  the  UAV  is  a  point  mass,  travelling  at  constant  velocity 
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Figure  2.1:  Pathmaker  solution  approach  [19] 


and  altitude,  with  a  constant  wind  vector  and  a  constant  target  velocity  along  its 
heading. 

Terning  also  showed  this  approach  is  rigorous  enough  to  handle  a  target  stopping 
and  changing  directions,  as  shown  in  Figure  2.2.  The  solid  black  line  on  the  left  is  the 
path  of  the  target,  the  blue  line  is  the  path  of  the  UAV.  The  target  moves  north  at 
10  m/s,  then  east  at  5  m/s,  south  at  20  m/s,  and  ending  west  at  15  m/s.  Also  in  this 
simulation  a  2  second  lag  is  incorporated.  As  long  as  there  is  enough  room  spacially 
to  track  the  target,  the  algorithm  can  accommodate  changes  in  the  target’s  path. 
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Figure  2.2:  Pathmaker  with  Target  Varying  Trajectory  [19] 

Figure  2.3  shows  there  are  instances  where  the  algorithm  has  trouble  intercept¬ 
ing  the  target.  Again,  the  target  path  is  black,  the  UAV  path  is  blue.  With  the 
combination  of  both  close  proximity  and  a  time  lag,  the  algorithm  can  display  oscilla¬ 
tory  behavior  and  fail  to  get  the  target  in  the  field  of  view  of  the  sensor.  This  would 
be  expected,  especially  with  a  UAV/target  speed  ratio  greater  than  2. 

Unfortunately  the  winter  weather  of  Ohio  prevented  Terning  from  rigorously 
flight  testing  this  concept  on  actual  hardware.  He  successfully  integrated  his  algorithm 
to  interface  with  Procerus  Technology’s  Kestrel  autopilot  system  (presented  in  Section 
3.1.2),  and  showed  in  simulation  that  the  concept  worked.  Although  not  truly  optimal, 
the  idea  appears  robust  enough  to  handle  communication  delays  (within  limits)  and 
quick  enough  to  reach  the  changing  target  conditions. 

This  result  shows  that  an  iterative  approach  to  changing  conditions  can  be 
robust  enough  for  real-time  applications.  This  approach  will  be  used  in  Chapter  111 
when  creating  the  cooperative  control  algorithm. 

2.2.5  Path  Planning -Rysdyk.  The  work  of  Rysdyk  [16]  presents  a  useful 

solution  to  keeping  a  target  in  sight  of  a  sensor  located  on  an  air  vehicle.  The  first 
part  uses  “helmsman  behavior”  control  scheme  to  keep  the  vehicle  on  the  desired 
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Figure  2.3:  Pathmaker  with  Large  Latency  [19] 

course.  The  second  part  develops  an  analytic  solution  to  the  correct  orbit  in  the 
presence  of  wind  to  keep  the  target  in  view.  “Helmsman  behavior”  determines  an 
appropriate  heading  to  get  back  on  the  desired  path  as  function  of  the  cross  distance 
(ys=perpendicular  distance  from  the  desired  path)  and  the  desired  heading  on  the 
path  (x.s=the  tangent  of  the  desired  path) (y  is  the  course/clock  angle  of  the  vehicle). 
The  equation  for  the  commanded  heading  becomes 

Xc(ys,Xs)  =  e(ys)  +Xs  (2.2) 

where  the  c  subscript  is  the  commanded  heading  and  a  is  a  function  that  behaves 
like  a  spring  (-kx)  but  saturates  at  the  values  [— Xicpt ,  Xicpt )]  (  Xicpt  denotes  intercept 
heading,  see  Figure  2.4).  The  author  chose  the  function 

e-ays/ 2  _  2 

v{ys)  =  Xicpt  e_aya/2  +  1  (2.3) 
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Figure  2.4:  Cross  Track  Function  [16] 


where  ys  is  the  crosstrack  distance  and  tilde  denotes  an  angle  with  respect  to  the 
desired  heading  (e.g.  \c  =  Xc  ~  Xs)-  If  you  choose  a  PID  control  scheme,  the  closed- 
loop  dynamics  becomes 


X  =  dy,Vasin(x)  +  kr(xc  -  x)  +  fc  /  (ic  -  x)d.T 


(2.4) 


where  Vg  is  the  ground  speed  and 

d  _  e~aVs/ 2 

ays  =  d^sa{ys)  =  ~aXlcp\e~^/ 2  +  l)2 
In  a  coordinated  turn,  bank  angle  is  related  to  turn  rate  so 


(2.5) 


tan(tl>)  =  -*-x 
9 


(2.6) 


This  can  be  used  to  provide  the  commanded  turn  rate  using  a  set  ground  speed 
(Vg)  and  is  a  fairly  simple  approach  to  maintaining  a  desired  path. 

The  above  analysis  assumes  constant  airspeed,  altitude,  and  coordinated  turns. 
The  reason  for  coordinated  turns  means  the  rate  of  change  of  the  heading  is  kine¬ 
matically  linked  to  the  bank  angle.  The  ideal  solution  makes  the  heading  exactly 
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tangent  to  the  desired  path  at  the  same  moment  the  cross-distance  becomes  zero. 
This  concept  will  be  used  to  model  closed-loop  UAV  behavior. 

The  concepts  in  both  previous  modelling  research  and  the  path  planning  back¬ 
ground  research  will  be  utilized  in  Chapter  III.  Many  ideas  like  the  above  (and  much 
more)  go  into  multi  UAV  control  algorithms.  The  following  “related  research”  is  a 
quick  survey  of  how  others  have  approached  the  problem  of  controlling  multiple  UAVs 
for  various  missions. 

2.3  Related  Research 

The  following  presents  4  different  approaches  to  cooperative  control  of  UAVs. 
The  scenarios  vary,  but  all  illustrate  the  complexity  of  this  area  and  also  illustrate 
the  pros  and  cons  of  optimality.  Suitability  of  the  algorithms  for  real-time  is  also 
presented. 

2.3.1  Broad  Area  Search-Beard.  Randal  Beard  [4]  explored  the  affect  of 
sub-optimal  solutions  on  both  performance  and  computational  time.  The  problem 
was  a  broad  area  search  with  randomly  generated  opportunities  and  threats.  The 
UAV  team  was  subject  to  the  constraints  of  vehicle  dynamics,  a  lower  distance  limit 
for  collision  avoidance,  and  an  upper  distance  limit  for  communications. 

One  example  used  3  alike  UAVs  with  a  front-looking  sensor  with  a  footprint  of 
width  w.  The  globally  optimal  solution  considered  all  possible  paths  which  maximized 
the  number  of  opportunities  observed  and  minimized  contact  with  the  threats.  In  this 
particular  scenario,  there  were  503  possible  paths.  This  brute  force  solution  sensed  10 
targets  and  took  522  seconds  to  solve.  By  constraining  the  optimal  feasible  vehicle 
paths  between  the  vehicle  path  on  the  right  and  the  vehicle  path  on  the  left,  10  targets 
were  observed  and  the  computational  time  decreased  to  13.3  seconds. 

The  sub-optimal  solutions  took  two  different  approaches.  The  first  approach 
lets  the  first  vehicle  take  the  “best  leader”  myopic  path  (not  considering  the  other 
vehicles).  The  second  vehicle  also  taking  the  next  myopic  path  only  constrained  by 
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the  first  path.  This  repeats  until  all  N  paths  are  generated  for  the  N  vehicles.  The 
second  sub-optimal  approach  takes  into  account  that  the  team  may  be  better  served  by 
each  individual  taking  sub-optimal  paths.  The  search  is  limited  to  “pairwise”  feasible 
paths.  With  the  “best  leader”  approach,  9  targets  were  sensed  and  the  algorithm 
found  a  solution  in  1.2  seconds.  The  author  did  not  run  the  second  sub-optimal  with 
the  3  UAV  example,  but  did  state  for  a  cost  function  J, 

Jbestleader  —  J pairwise  A  3 optimal  (2.7) 

This  article  is  a  good  cost-benefit  analysis  between  optimal  and  sub-optimal 
solutions,  and  similar  results  are  expected  for  the  heuristic  approach  developed  herein. 

2.3.2  Task  Oriented  Cooperative  Control-  Richards /Bellingham.  Arthur 

Richards  and  John  Bellingham  [15]  of  MIT  added  a  couple  more  layers  of  complexity 
to  the  problem.  The  problem  was  constructed  to  accomplish  a  mission  in  minimum 
time  (the  cost  function).  The  “mission”  became  visiting  all  required  waypoints  and 
also  adhering  to  any  timing  constraints  (e.g.  must  visit  A  five  minutes  after  D).  Not 
all  UAVs  were  alike,  and  only  certain  vehicles  could  visit  certain  waypoints  depend¬ 
ing  on  their  capabilities.  Superimposed  on  this  was  “no  fly  zones”  -rectangular  areas 
where  trajectories  could  not  penetrate. 

The  intent  of  this  method  was  to  find  a  globally  optimal  solution.  Realizing  the 
computational  intensive  nature  of  this  solution,  the  authors  suggest  using  this  as  a 
benchmark  on  which  heuristic  methods  can  be  compared,  and  is  not  well  suited  for 
real-time  applications.  To  illustrate  the  optimal/approximate  differences,  the  solution 
is  compared  to  a  less  intensive  method  that  estimates  the  trajectory  planning  portion 
of  the  problem. 

The  UAV  model  was  simplified  to  two  dimensions.  The  vehicle  itself  was  mod¬ 
elled  as  a  point  mass.  The  x-y  position,  x-y  velocity,  and  yaw  rate  comprised  the 
whole  state.  Vehicle  dynamics  were  linearly  modelled.  The  velocity  is  constrained  to 
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a  maximum  magnitude,  and  the  turn  rate  is  also  constrained  by  placing  a  limit  on 
the  lateral  force  magnitude.  The  optimization  cost  function  minimizes  not  only  time 
but  also  weights  control  effort  of  each  UAV. 

The  approximate  method  simplifies  the  solution  search.  First,  assignments  that 
place  a  large  number  of  waypoints  (tasks)  onto  a  single  vehicle  are  eliminated.  To 
find  approximate  completion  times,  a  straight  line  approximation  method  is  used. 
Once  costs  are  calculated  and  an  approximate  minimum  cost  is  found,  then  detailed 
trajectory  analysis  is  performed  that  accounts  for  dynamics  and  collision  avoidance. 

Optimal  and  sub-optimal  results  were  compared  with  a  scenario  involving  three 
UAVs,  2  no  fly  zones,  and  4  waypoints.  Using  CPLEX1  optimization  software  with  a 
1GHz  PC  with  256MB  of  RAM,  the  computation  took  over  eight  minutes.  With  the 
same  PC,  CPLEX  and  MATLAB  software,  the  approximate  method  found  a  solution 
in  under  5  seconds.  In  this  case,  the  approximate  method  found  the  globally  optimal 
solution. 

This  is  a  good  illustration  that  when  approximations  are  used  in  the  appropriate 
locations,  solutions  can  be  nearly  optimal  (maybe  even  optimal)  at  a  fraction  of  the 
computer  time.  This  research  is  a  variation  of  Beard’s  approach,  but  requires  that 
more  than  1  UAV  is  needed  to  accomplish  a  task.  When  compared  to  optimal,  the 
algorithm  presented  in  this  research  will  also  be  comparable  to  the  “best  leader” 
approach,  since  it  also  takes  a  myopic  approach  based  on  1  UAV  (See  Chapter  3). 

2.3.3  Cooperative  Control  Simulation- Rasmussen/ Mitchell.  To  deploy  a  us¬ 
able  system  that  utilizes  UAVs  in  a  cooperative  manner,  many  factors  besides  optimal 
path  planning  affect  the  performance.  Whether  the  mission  is  a  broad  area  search, 
surveillance/reconnaissance,  or  search  and  destroy,  the  necessary  tasks  may  require 
UAVs  to  detect,  assess,  engage,  assess  again,  and  possibly  engage  again.  Any  one  of 
these  tasks  are  not  trivial  for  an  automated/ unmanned  system  (let  alone  a  manned 

1A  software  package  sold  by  ILOG(IBM)  for  business  decision  making,  efficient  resource  utiliza¬ 
tion,  and  scheduling/planning.  For  more  info:  www.ilog.com/products/cplex 
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system).  AFRL’s  Air  Vehicles  directorate  created  a  simulation  environment,  Mul- 
tiUAV2  [14],  that  attempts  to  provide  an  assessment  tool  for  cooperative  approaches, 
but  is  not  intended  for  real-time  use.  The  thought  that  went  into  this  product  illus¬ 
trates  the  complexity  of  the  problem,  and  highlights  items  to  consider  as  an  engineer 
or  researcher. 

The  state  of  all  targets  is  passed  to  all  vehicles.  Each  target  state  includes 
position,  and  whether  the  team  has  detected,  classified,  attacked,  and  assessed  the 
target.  Messages  are  passed  to  all  as  a  target  changes  state.  Each  vehicle  computes 
the  cost  to  accomplish  the  remaining  tasks,  and  transmits  to  all.  The  cooperative 
control  algorithm  resides  on  all  vehicles,  the  cost  is  assessed  and  assignments  are 
given. 

At  the  lowest  level  of  logic  is  vehicle  dynamics  and  path  planning.  All  vehicles 
are  modelled  with  a  6  degree  of  freedom  aircraft  model,  inner  control  loops  that 
govern  the  desired  velocities  and  attitude,  and  outer  control  loops  that  govern  desired 
altitude,  heading,  and  waypoint  tracking.  One  step  higher  in  the  hierarchy  is  the 
path.  Unless  specified  otherwise,  all  paths  calculated  are  optimal.  The  task  at  the 
end  of  the  path  could  take  many  forms:  a  post-attack  assessment,  an  attack,  a  second 
look,  or  possibly  an  anomaly  that  needs  classification. 

The  detection  portion  combines  probability  and  directional  dependence.  This 
accounts  for  the  fact  that  even  if  a  target  is  in  the  field  of  view,  there  is  no  guarantee 
that  it  will  be  detected.  The  probability  of  identification  varies  depending  on  which 
direction  the  sensor  views  the  target.  Using  trigonometric  functions,  the  probability 
maximizes  at  defined  directions  and  decreases  as  the  view  angle  deviates  from  the  ideal 
angle.  For  example,  in  the  default  scenario,  for  a  single  pass,  the  maximum  certainty 
is  0.8  and  you  need  a  certainty  of  0.9  before  you  can  attack.  So  the  simulation  forces 
multiple  passes  before  the  target  can  be  positively  identified. 

At  the  highest  level  are  the  Cooperation  algorithms.  Each  vehicle  has  these  al¬ 
gorithms  on-board,  and  all  vehicles  are  “in  sync”  with  their  assigned  tasks.  Given  the 
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state  of  the  system  (including  all  UAV  states  and  all  target  states),  tasks/assignments 
are  assigned  to  each  vehicle.  The  different  algorithms  to  distribute  the  workload  vary 
greatly,  but  fall  into  two  general  categories:  Single  Assignment  Tour  and  Multiple 
Assignment  Tour.  As  the  name  implies,  Single  Assignment  Tour  hands  the  vehicle 
only  one  task  (e.g.  search,  or  classify)  at  a  time.  It  keeps  complexity  down,  but  can 
be  very  inefficient.  The  Multiple  Assignment  Tour  accounts  for  the  next  step  in  the 
process  and  improves  efficiency,  but  can  lead  to  a  “combinatorial  explosion.”  Many 
algorithms  can  be  implemented,  both  optimal  and  suboptimal.  Brief  explanations  for 
each  approach  are  included  in  [14]. 

This  tool  forces  any  designer  of  a  Cooperative  Control  scheme  to  consider  the 
necessary  attributes  needed  to  create  a  useful  product  to  the  warfighter.  The  de¬ 
sign  must  include  much  more  than  path  planning  and  should  consider  probability 
of  detection,  complexity  vs  efficiency,  and  inter  UAV  communication.  A  real-time 
autonomous  application  of  this  extends  well  beyond  the  scope  of  this  research,  and 
provides  a  road  to  follow  for  future  real-time  research. 

2.4  Summary 

As  this  chapter  illustrates,  using  UAVs  in  a  cooperative  manner  pulls  knowledge 
from  many  areas  and  often  becomes  a  very  complex  problem.  To  aid  in  tackling  this 
task,  this  research  pulls  background  information  from  past  work  in  aircraft  modelling 
and  path  planning.  The  works  of  Jodeh  and  Kingston  will  aid  in  creating  a  repre¬ 
sentative  model.  Also  the  works  of  Zollars  and  Terning  establish  a  base  to  create 
a  path  planning  approach.  The  related  research  section  illustrates  both  the  varied 
application  of  multi-UAV  schemes  and  also  the  trade-offs  of  optimal  vs  sub-optimal 
solutions.  This  provides  a  good  base  to  create  a  real-time  multi  vehicle  algorithm. 
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III.  Equipment  and  Analysis 


3. 1  Introduction 

This  chapter  lays  the  foundation  for  modelling  the  BATCAM  and  creating  the 
cooperative  control  algorithm.  It  begins  with  a  description  of  the  BATCAM  and  the 
equipment  associated  with  flight  testing.  Pulling  from  research  presented  in  Chapter 
II,  a  closed-loop  model  is  created  for  the  BATCAM.  Next  the  algorithm  is  developed 
to  simultaneously  approach  and  maintain  a  surveillance  orbit.  The  chapter  ends  with 
a  brief  overview  of  the  cooperative  control  algorithm  software  that  interacts  with  each 
UAV. 

3. 2  Equipment 

3.2.1  BATCAM.  Applied  Research  Associates  (ARA)  in  conjunction  with 
AFRL  developed  the  BATCAM.  The  idea  came  from  the  need  to  develop  a  tactical  air 
surveillance/reconnaissance  tool  with  easy  portability  and  very  low  logistical  needs. 
The  primary  customer  for  this  product  are  special  operations  forces.  The  BATCAM 
provides  real  time  situational  awareness  and  targeting  information  [3]  and  has  many 
features  well  suited  for  field  use  (illustrated  in  Figure  3.1).  Both  its  light  (0.85  lb) 


Figure  3.1:  BATCAM  UAV  [3] 

carbon  composite  airframe  and  flexible  composite/fabric  wings  makes  the  body  re¬ 
silient  to  damage.  Launching  is  done  by  hand  (like  a  paper  airplane)  and  is  powered 
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Table  3.1:  BATCAM  Camera  FOV  Angles 


Angle 

Front  Camera 

Side  Camera 

Depression  Angle 

49° 

39° 

horizontal  FOV 

O 

00 

O 

00 

vertical  FOV 

o 

O 

o 

O 

by  a  quiet  electric  motor.  The  batteries  can  be  recharged  with  a  car  cigarette  lighter 
in  less  than  an  hour.  Assembly  of  the  wing  airframe  takes  less  than  a  minute,  and 
removing  a  spent  battery  with  a  recharged  battery  also  takes  about  a  minute. 

The  surveillance  capability  consists  of  two  small  cameras  hard  mounted  to  the 
bottom  of  the  fuselage.  The  first  camera  is  a  forward  looking  camera  ,  the  second 
camera  is  a  side  looking  camera  facing  left  (from  the  pilot’s  perspective).  Neither 
camera  is  gimballed,  so  the  attitude  of  the  BATCAM  dictates  the  FOV  of  the  sensor. 
Table  3.1  provides  the  specifics  on  the  FOVs.  The  side  camera  will  be  used  for  the 
algorithm  as  the  UAV  orbits  the  target  counter-clockwise. 

3.2.2  Kestrel  Autopilot.  Each  BATCAM  used  for  flight  test  contained  Pro¬ 
cerus  Technology’s  Kestrel  Autopilot.  This  small  electronic  device  (see  Figure  3.2) 
provides  “autonomous  flight  control,  Global  Positioning  System  (GPS)  waypoint  nav¬ 
igation,  autonomous  take-off  and  landing.”  [1]  The  sensor  suite  contains  3-axis  rate 
gyros,  accelerometers,  differential  and  absolute  air  pressure  sensors  for  attitude  esti¬ 
mation  and  altitude/airspeed  estimation.  A  GPS  receiver  provides  positional  infor¬ 
mation  to  the  flight  computer.  Extra  serial  ports  provide  the  ability  to  dynamically 
communicate  and  execute  commands  to  cameras  or  payload  devices.  The  autopilot 
communicates  with  the  ground  station  at  900  MHz  with  specialized  TCP/IP  packets. 
In  addition  to  command  packets,  telemetry  packets  provide  the  aircraft  state  infor¬ 
mation  at  rates  up  to  25  Hz.  The  autopilot  is  just  one  component  of  the  complete 
system. 

The  complete  system  requires  a  laptop  for  human  interface,  software  for  the 
laptop,  a  USB  communication  box  to  “talk”  to  the  autopilot,  the  autopilot  itself,  and 
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Figure  3.2:  Kestrel  Autopilot  [1] 

a  Radio  Controlled  type  model  airplane  retrofitted  with  a  pitot  tube.  The  software, 
named  Virtual  Cockpit,  provides  the  user  the  ability  to  command  a  loiter,  a  set  of 
waypoints,  an  airspeed,  an  attitude,  and  set  up  fail  safes  in  the  event  of  problems.  VC 
can  also  toggle  between  manual  and  autonomous  modes,  monitor  progress  with  de¬ 
tailed  displays,  and  adjust  many  parameters  in-flight.  The  software  can  communicate 
to  multiple  UAVs  simultaneously  via  the  communications  box.  The  Virtual  Cockpit 
(VC)  software  will  be  the  primary  means  to  command  the  autopilots  in  the  air. 

The  software  development  kit  for  Virtual  Cockpit  creates  an  interface  to  create 
custom  applications.  Since  the  communication  scheme  is  TCP/IP  based,  all  informa¬ 
tion  going  to  and  from  the  UAVs  is  packet  based.  Many  packet  types  are  available, 
some  that  command  the  autopilot,  and  some  contain  UAV  telemetry.  Using  the  stan¬ 
dard  loop-back  network  capability,  the  custom  application  can  run  alongside  Virtual 
Cockpit,  receive  all  incoming  packets,  and  pass  command  packets.  This  is  the  ap¬ 
proach  used  to  demonstrate  the  control  scheme. 
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3.2.3  Flight  Test  Setup.  The  equipment  used  for  flight  testing  is  a  self 
contained  mobile  trailer  with  all  the  necessary  items  for  both  flying  and  repair.  The 
trailer  has  both  an  operations  area  and  a  maintenance/repair  area  (see  Figure  3.3). 

The  operations  area  has  all  the  necessary  capability  to  both  command  and  con¬ 
trol  the  UAVs  and  display  video.  The  trailer  provides  the  operators  and  electronics 
protection  from  the  elements.  The  laptop  is  used  to  command  and  control  the  UAVs 
and  the  video  equipment  provide  the  surveillance  capability.  The  laptop  passes  com¬ 
mands  and  receives  telemetry  from  the  comm  box  using  a  serial  to  USB  connection. 
This  comm  box  is  also  made  by  Procerus  Technologies  (who  makes  the  autopilot)  and 
is  a  necessary  component  to  communicate  with  the  autopilot.  The  900  MHz  antenna 
is  connected  to  the  comm  box.  Video  is  received  via  two  2.1  GHz  antennas.  The 
antennas  are  configured  to  maximize  coverage.  Each  signal  is  split  into  4  and  passed 
to  4  receivers,  providing  the  capability  to  receive  4  separate  feeds  from  4  UAVs.  Each 
receiver  has  two  feeds,  and  chooses  the  one  with  the  best  quality  for  the  output.  The 
video  switch  takes  the  4  feeds  and  provides  the  flat  screen  TV  a  customizable  display. 
Possible  displays  include  all  4  feeds  on  1  screen,  two  chosen  feeds,  or  just  1  feed  from 
any  available  signal.  Both  the  command,  control  and  surveillance  capabilities  provide 
a  good  platform  for  UAV  research. 


Top  View 


Figure  3.3:  Flight  Testing  Setup 
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The  maintenance/repair  area  helps  sustain  the  UAVs.  The  toolbox  contains  a 
variety  of  basic  tools  (wrenches,  screwdrivers,  socket  sets,  glues)  and  specialized  equip¬ 
ment  for  Radio  Controlled  (RC)  airplanes  (e.g.  refuelling  pumps,  engine  starters). 
The  workbench  provides  an  area  to  work  and  also  has  a  strip  of  outlets  for  battery 
rechargers  and  power  tools.  Also  on  hand  are  repair  materials  like  balsa  wood,  tubing, 
plastic,  extra  propellers,  etc. 

Flight  Tests  are  conducted  at  Camp  Atterbury,  Indiana.  Since  the  Federal 
Aviation  Administration  (FAA)  requires  that  autonomous  vehicles  be  only  flown  in 
Restricted  Airspace,  Camp  Atterbury  is  the  nearest  facility  to  AFIT  with  this  desig¬ 
nation.  The  UAVs  are  flown  at  the  airfield,  and  each  flight  is  coordinated  with  the 
airfield. 

Operations  are  conducted  in  typical  Air  Force  style.  Test  objectives  are  created 
for  each  flight  test  and  test  cards  are  written.  Each  card  contains  detailed  procedures 
necessary  to  accomplish  each  test.  Before  testing,  the  whole  test  is  presented  to  a 
Technical  Review  Board/Safety  Review  Board  (TRB/SRB)  for  approval.  After  incor¬ 
poration  of  comments  and  approval  from  the  presiding  officer,  the  test  is  conducted. 

Personnel  are  assigned  specific  roles  to  accomplish  the  test.  The  Test  Conduc¬ 
tor  oversees  the  test  operations.  The  Pilot  is  a  certified  UAV  operator  that  keeps  the 
vehicles  in  sight  at  all  times  and  can  manually  fly  them  if  necessary.  The  Ground 
Operator  sits  at  the  laptop  and  controls  the  UAVs  using  the  VC  software.  The  Data 
Recorder  documents  all  information  for  each  test  card  and  also  any  anomalies  noted. 
The  Launch/Recovery  personnel  accomplish  launch  and  recovery.  The  Safety  Officer 
supervises  the  test  and  provides  input  to  minimize  injury  and  equipment  damage. 
During  flight,  the  Ground  Operator  and  Pilot  are  in  constant  communication  so  that 
the  Pilot  understands  what  the  UAV  is  suppose  to  be  doing.  The  Pilot  also  pro¬ 
vides  feedback  to  the  Operator  to  improve  performance  and  prevent  any  undesired 
situations. 
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3.3  Airframe- Properties 

The  BATCAMs  came  to  AFIT  “as  is.”  To  make  up  for  the  lack  of  technical 
drawings  and  sparse  documentation,  some  basic  measurements  were  needed  to  begin 
the  modelling  process.  Physical  dimensions  and  mass  were  combined  with  Kingston’s 
[9]  closed-loop  model  to  create  a  general  airframe  description. 

3.3.1  Physical  Dimensions.  Airframes  are  typically  geometrically  complex, 
and  the  BATCAM  is  no  exception.  The  desired  approach  needed  to  be  simple,  effi¬ 
cient,  but  also  reasonably  accurate  and  complete.  For  the  purposes  of  this  research, 
highly  accurate  measurements  were  not  necessary.  For  dimensions  not  easily  measured 
with  a  tape  measure,  a  photographic  technique  was  used.  Digital  pictures  are  taken 
of  the  airframe  from  each  of  the  major  axes  (side,  front,  rear,  top),  with  a  reference 
measurement  visible  in  the  picture.  To  minimize  parallax  and  perspective  errors,  the 
reference  measurement  must  be  placed  at  the  same  distance  as  the  item  of  interest, 
the  focal  axis  must  be  orthogonal  to  the  image  plane,  and  the  distance  between  item 
and  camera  must  be  sufficient.  This  technique  is  regularly  used  in  repair  situations 
where  the  engineer  is  physically  separated  from  the  damaged  part  but  must  design  a 
fix. 

Figure  3.4  illustrates  this  technique.  Using  software  comparable  to  GNU  Image 
Manipulation  Program  (GIMP),  pixel  distances  can  be  measured.  In  this  example, 
the  desired  measurement  is  the  overall  length  of  the  BATCAM.  Three  pixel  mea¬ 
surements  are  required:  the  pixel  distance  of  1  inch(ppi),  the  pixel  distance  of  the 
BATCAM(dpja;),  and  the  angle  of  the  bottom  of  the  UAV  (9 bat)  with  respect  to  the 
picture.  The  approximate  length  (1)  becomes 

l  =  ~  dp‘‘  ,  (3.1) 

ppi  COS (Ubot) 

In  this  case,  the  horizontal  pixel  distance  is  1965  pixels,  ppi  is  98  pixels  per  inch,  and 
the  bottom  of  the  BATCAM  is  tilted  1.54  degrees,  yielding  an  overall  length  of  20.06 
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Figure  3.4:  Length  Measurements  using  Digital  Pictures 


inches.  From  the  BATCAM  users  manual  [3],  the  stated  length  is  20  inches.  When 
done  properly,  the  margin  of  error  for  this  technique  is  estimated  to  be  between  .2 
and  1%. 

3.3.2  Mass,  Aerodynamics,  Results.  This  section  presents  the  airframe  prop¬ 
erties  required  for  the  model  in  section  3.4.1.  It  includes  a  set  of  physical  dimensions 
and  areas  needed  for  lift  calculations  and  mass  properties. 

The  mass  is  determined  by  simply  placing  the  BATCAM  on  a  scale.  From  basic 
physics,  the  weight  is  divided  by  the  gravitational  acceleration  to  find  mass.  The  scale 
used  was  accurate  to  the  nearest  gram. 

In  Section  3.4.1,  the  only  needed  aerodynamic  force  for  the  model  is  lift.  This 
research  will  assume  a  constant  coefficient  (Cl)  of  lift  at  0  degrees  angle  of  attack 
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Table  3.2:  Dimensions  and  Mass  Properties 


Item 

Value 

Wing  Reference  Area  (S) 

103.7  in2 

Mass  (m) 

0.425  kg 

(AoA)  centered  around  the  level  steady  flight  conditions.  From  flight  testing,  level 
steady  flight  occurred  at 


Vao  =  11.75  m/s 
Throttle  =  50% 


(3.2) 


DeLuca  [5]  conducted  an  extensive  aerodynamic  analysis  of  the  BATCAM  airframe 
at  airspeeds  of  10,  20,  30,  and  50  mph.  Based  on  the  steady  flight  conditions,  the  20 
mph  data  for  a  flexible  was  used  to  determine  Cl-  Using  basic  fluid  mechanics  [12], 
lift  (L)  for  the  BATCAM  is 


L 


Cl  =  0.7 


(3.3) 


where  Va  is  airspeed,  p  is  the  air  density,  and  S  is  the  wing  reference  area.  All  other 
aerodynamic  forces  are  not  omitted,  but  are  rolled  into  the  first  order  constants  in 
the  model.  Section  3.4.1  will  explain  this  further.  Table  3.2  summarizes  the  needed 
values  for  the  model. 


3-4  Simplified  Closed  Loop  model 

This  section  creates  a  model  for  a  single  UAV  that  describes  the  closed  loop 
system  created  by  the  UAV  and  autopilot.  As  stated  in  the  assumptions,  the  response 
of  some  of  the  state  is  assumed  to  mimic  a  first-order  ordinary  differential  equation 
(ODE),  discarding  higher  order  terms.  Other  variables  in  the  state  retain  the  non¬ 
linear  properties. 


3-4-1  Base  Model.  The  simplified  model  is  based  upon  Equation  2.1.  To 
account  for  dynamics  in  the  vertical  direction,  pitch  ( 6 )  and  altitude  (h)  become  part 
of  the  state,  and  Equation  2.1  is  changed  to 


Pn  —  Va  cos  -0  cos  6  +  Wn 
Pe  =  Va  sin  0  cos  9  +  we 


Va 

h 

h 

0 

9 

0 


kv(Vac  -  Va) 


Va  sin  9 

L(Va)  cos  0  cos  9 
m 

£^(0C  -  0)  +  W'p 

ke(9c-9) 


(3.4) 


For  this  set  to  be  reasonably  accurate,  roll  (0)  and  pitch  (9)  must  remain  “small”. 
The  proportional  gains  k*  will  be  empirically  calculated  from  flight  test  results.  Note 
that  changing  the  heading  depends  on  banking  (roll)  the  UAV  (i.e.  a  coordinated 
turn).  The  commanded  pitch  9C  will  be  used  to  reach  the  desired  altitude.  The  wind 
disturbance  ( w <p)  in  the  roll  rate  equation  captures  the  impact  of  wind  disturbances 
on  the  side  camera  of  the  BATCAM.  Many  of  the  aerodynamic  forces  (except  L)  are 
included  in  the  k*  values.  By  using  flight  test  data  centered  around  the  nominal  cruise 
airspeed,  the  extracted  k  values  should  reasonably  reflect  BATCAM/autopilot  system 
behavior. 
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The  wind  will  be  modelled  using  the  Dryden  wind  turbulence  model  [10].  This 
representation  adds  turbulence  to  velocity  spectra  by  passing  band  limited  white  noise 
through  appropriate  filters.  These  filters  are  mathematically  described  in  MIL-HDBK 
1797  [10].  The  roll  disturbance  will  be  a  scaled  quantity  of  the  wind  magnitude.  Con¬ 
veniently,  MATLAB’s  Simulink  program  has  a  prepackaged  block  that  produces  this 
wind  model  in  the  Aerospace  blockset.  The  inputs  to  this  block  are  wind  speed  and 
direction  at  6  meter  height,  a  probability  scale  for  light/medium/heavy  turbulence, 
scale  height,  UAV’s  airspeed,  UAV’s  altitude,  and  the  output  is  the  North/East /Down 
components  of  wind.  Only  the  North  and  East  components  are  used  in  this  research. 
The  roll  disturbance  is  a  scaled  value  of  the  magnitude  of  the  wind  (w^  =  \  W  | ) . 
Figure  3.5  is  a  typical  wind  profile  for  a  1  m/s  input. 


Figure  3.5:  Dryden  Wind  Model-1  m/s  Input  Wind  Profile 

A  full  6  DOF  model  based  on  measured  data  was  initially  attempted,  but  not 
yet  completed.  Appendix  A  contains  this  more  rigorous  model  of  the  BATCAM  (and 
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included  MATLAB  files)  that  includes  linearized  aerodynamic  coefficients  based  on 
DeLuca’s  work  [5]. 


3-4-2  Closed  Loop  Model  Details.  Both  commanded  heading  (yc)  and  com¬ 
manded  altitude  (hc)  are  key  parameters  in  controlling  the  UAV,  but  neither  of  these 
appear  in  the  model  above.  Changes  in  heading  depend  on  roll  angle,  and  changes 
in  altitude  depend  on  airspeed  and  pitch  angle.  This  section  provides  the  details  to 
bridge  roll  to  heading,  and  altitude  to  pitch. 

Cross  distance  (ys)  is  the  perpendicular  distance  of  the  UAV  from  the  intended 
path  (See  Fig  2.4  or  Fig  3.6).  This  is  an  important  quantity  for  two  reasons:  1  -ys 
will  be  used  as  a  measure  of  performance  for  the  algorithm;  2 -ys  is  used  to  determine 
the  commanded  heading  yc.  Since  UAV  placement  is  one  of  the  key  factors  affecting 
sensor  performance,  Chapter  IV  will  use  ys  to  quantify  algorithm  performance.  ys 
is  also  used  to  provide  feedback  and  get  the  UAV  back  on  the  desired  path.  The 
following  analysis  will  find  the  cross  distance  using  basic  vector  math  and  then  use 
Rysdyk’s  [16]  approach  presented  in  Section  2.2.5  to  find  the  commanded  heading  to 
get  the  UAV  back  on  the  desired  path.  The  error  between  commanded  heading  and 
UAV  heading  will  feed  into  the  roll  feedback  equation  to  cause  a  turn  rate  ip  (see 
Equation  3.2). 

Figure  2.4  illustrates  that  the  commanded  heading  depends  on  both  the  heading 
of  the  desired  path(ys)  and  the  cross  distance  from  the  desired  path(r/s).  The  first 
step  to  find  x°  is  t°  find  Vs •  Figure  3.6  illustrates  the  vectors  necessary  to  fold  ys. 
The  waypoint  location  is  at  (Pe2,  Pn2)—P2,  the  desired  straight  line  path  starts  at  the 
origin  to  (Pe2,  Pm),  and  the  current  UAV  location  is  at  (Pe,  Pn)—P-  The  desired 
behavior  is  for  the  UAV  to  go  to  the  waypoint  along  the  desired  path.  From  the 
definition  of  the  cross  product,  the  crosstrack  distance  (ys)  is 


Vs 


P2  X  P 
\p2\ 


(3.5) 
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Realize  that  ys  is  positive  when  9\  <  02  and  ys  is  negative  when  9i  >  92.  To  drive  ys 
to  zero,  the  commanded  heading  equation  (Eqn  2.2)  becomes 


X  =  ^ 


7 T 


\P2  X  P | 
\Po\ 


Xs  =  --  arctan 


+  X 

PN2 

Pe2 


(3.6) 


The  a  function  was  defined  in  Equation  2.3.  Also  note  that  the  second  equation 
converts  an  angle  using  the  mathematical  convention  to  a  clock  angle/heading  angle. 
The  commanded  roll  angle  becomes 


kx(a(Vs)  +  arctan  -  x)  lxc  ~  xl  <  Xr: 


(3.7) 


|XC  -  Xl  >  Xr, 


When  the  heading  is  “close”  to  the  commanded  heading,  then  a  proportional  gain 
will  be  used.  When  the  heading  difference  is  large  enough,  the  commanded  roll  angle 
(which  controls  turn  rate)  is  limited  to  a  set  value  (< f>max )•  This  is  to  keep  the  UAV 
stable.  The  value  Xmax  will  determine  this  threshold. 

The  commanded  pitch  depends  on  the  error  between  commanded  altitude  hc 
and  UAV  altitude  h.  Bridging  altitude  to  pitch  will  also  use  a  proportional  gain.  The 
commanded  pitch  will  be 


6C  = 


kh(hc  —  h )  | hc  —  h\  <  hn 

9 max  \h  h\  >  hr 


(3.8) 


Just  like  roll,  commanded  pitch  will  reach  a  maximum  value  when  the  altitude  differ¬ 
ence  exceeds  a  certain  value  ( hmax ).  Once  again,  this  is  to  maintain  positive  control 
of  the  UAV. 

Up  to  this  point,  the  analysis  has  only  described  the  BATCAM/autopilot  model, 
and  has  not  started  the  cooperative  control  algorithm.  The  inputs  for  this  model  are 


32 


Figure  3.6:  Illustration  of  Calculating  Crosstrack  Distance 

the  UAV  state,  waypoint  location  (pjv2,  PE2),  the  desired  path  heading  Xs,  com¬ 
manded  altitude  hc,  and  commanded  airspeed  I/c.  The  cooperative  control  algorithm 
must  provide  these  inputs  (except  the  UAV  state).  Section  3.6  will  describe  how  the 
algorithm  will  provide  these  inputs  based  on  the  state  of  all  UAVs.  One  last  section  of 
background  analysis  remains  for  describing  the  individual  UAV  before  the  cooperative 
control  algorithm,  the  Sensor  Footprint. 
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3.5  Sensor  Footprint  Analysis 

This  section  creates  a  projection  of  the  FOV  onto  the  ground  given  the  UAV’s 
attitude.  As  stated  in  Chapter  I,  this  analysis  assumes  that  the  ground  is  flat  and 
that  the  primary  sensor  is  the  side  camera.  This  analysis  provides  the  means  to 
determine  if  a  known  target  is  visible,  and  also  provides  a  way  to  determine  the 
surveillance  performance  of  the  UAV.  If  the  target  lies  within  the  FOV  quadrilateral, 
then  the  target  is  “visible”  to  the  camera.  The  measure  of  performance  will  find  the 
percentage  of  time  while  orbiting  that  the  target  is  visible. 

The  approach  begins  by  creating  vectors  in  the  body  frame  and  rotating  them 
into  the  inertial  frame  via  a  directional  cosine  matrix.  Then  a  parametric  line  will 
be  drawn  from  the  UAV  position  to  the  ground  (z=0)  in  the  direction  of  the  rotated 
vector. 


Front  View 


Top  View 


Figure  3.7:  BATCAM’s  Field  of  View 

3.5.1  FOV  from  body  frame  to  inertial  frame.  Please  refer  to  Figure  3.7. 
The  BATCAM’s  FOV  for  the  left  facing  side  camera  is  square  with  the  angles  shown. 
Each  corner  of  the  FOV  must  be  turned  into  a  vector(e)  in  M3.  If  is  the  depression 
angle,  FOVh  is  the  horizontal  FOV,  FOVy  is  the  vertical  FOV,  then  the  vectors  for 
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all  4  corners  of  the  FOV  is 


(  tan  F°v»  \ 

tan  F°p  ) 

-1 

£2  = 

-1 

ytan  (7 d  +  FC ^  )  J 

\tan  (7 d  FC, y  ) ) 

tan  77  ) 

tan  77  \ 

e\  = 

-1 

et  = 

-1 

ytan  (7d  +  fo2h  ) ) 

\tan  (7 d  F°^  )  / 

(3.9) 


These  vectors  in  the  body  frame  use  the  standard  body  frame  convention  where  the 
origin  is  at  the  center  of  gravity,  the  +x-axis  extends  out  the  nose  of  the  UAV,  the  +y 
axis  extends  out  the  right  wing  (from  the  pilot’s  perspective),  and  +z-axis  extends 
down.  The  b  denotes  body  frame.  Note  the  negative  values  in  the  y  component  to 
indicate  a  left  facing  camera.  To  rotate  any  of  these  body  vectors  into  the  geographic 
North/ East /Down  (NED)  frame  a  directional  cosine  matrix  (DCM)(C'^1  [18])  is  used 
that  utilizes  the  roll(</>),  pitch(0)  and  yaw(r/)  angles  (Euler  Angles).  The  vector  in 
the  NED  frame  is 

^1,2, 3, 4  =  Cb  ^1,2, 3, 4  (3.10) 

where  (using  the  abbreviations  c=cos,  s=sin) 


CT 


cOcip  —ceps'll)  +  sepsOc'ip  seps'ip  +  cepsOcip 
cos'll  expe^  +  sepsOsip  —sepc'ip  +  cepsOsip 
—s6  sepcO  cepcO 


(3.11) 


The  n  denotes  the  geographic  NED  frame.  Now  the  vectors  are  ready  to  extend  to 
the  ground. 


3.5.2  Creating  the  Footprint.  With  the  UAV  at  position  ( Pn,Pe,Pd )  and 
attitude  ( ep ,  0,  ip),  lines  that  pass  through  the  position  with  slopes  e  need  to  be  created. 
Parametric  representation  is  a  convenient  representation  for  this  purpose.  Using  the 
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parameter  (t),  not  to  be  confused  with  time,  the  vector  representation  of  any  point 
on  the  line  is 


H 

yU 

—  0,2, 3, 4^  T 

Pe 

V) 

[pd) 

(3.12) 


These  lines  lie  on  each  corner  of  the  FOV  and  extend  from  the  UAV  position  to  the 
ground.  There  will  be  a  total  of  4  lines  for  the  square  FOV. 


Figure  3.8:  Sample  Field  of  View  Projected  onto  the  Ground 

It  is  important  to  remember  that  in  the  NED  convention,  a  height  above  the 
ground  shows  up  as  a  negative  number  since  +z  is  down,  therefore  altitudes  should  be 
negative.  Taking  the  zn  component  and  solving  for  t  when  zn= 0  provides  the  point  on 
the  line  where  the  FOV  intersects  the  ground.  Each  of  the  four  e  vectors  will  create 
a  point  on  the  ground  (z=0  plane).  By  drawing  a  line  between  each  of  these  points, 
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a  quadrilateral  is  formed.  Figure  3.8  illustrates  this.  If  the  target  resides  within  this 
quadrilateral,  the  target  is  in  the  FOV.  There  are  instances  where  the  quadrilateral 
cannot  be  drawn.  If  the  attitude  of  the  UAV  causes  the  FOV  to  be  above  the  horizon, 
then  a  solution  does  not  exist.  The  MATLAB  routine  written  to  create  the  FOV  on 
the  ground  will  check  for  this  condition. 
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3.6  Cooperative  Control  Algorithm 

This  algorithm  was  written  to  command  multiple  UAVs  to  orbit  a  single  tar¬ 
get.  The  goals  for  the  algorithm  are  to  have  all  UAVs  enter  the  surveillance  orbit 
simultaneously,  then  maintain  equal  angular  spacing  during  the  orbit.  This  algorithm 
assumes  the  UAV  is  using  a  side  camera  and  views  the  target  while  maintaining  a 
CCW  circular  path  above  the  target. 

3.6.1  Overview.  This  CC  approach  creates  Dubins  Paths  [6]  for  each  UAV 
from  their  initial  position  and  heading  to  a  final  position  and  heading  that  enters  a 
surveillance  orbit  all  tangent  to  the  orbit  and  all  in  the  same  direction  (CCW).  To 
coordinate  entering  the  surveillance  orbit  in  an  equally  spaced  manner,  each  path  for 
each  UAV  is  generated  such  that  all  have  equal  lengths,  and  are  commanded  such 
that  the  ground  speed  is  equal  for  all.  To  aid  each  UAV  so  that  they  all  enter  the 
surveillance  orbit  at  the  same  time,  velocity  feedback  is  incorporated  as  they  traverse 
the  path.  As  all  UAVs  enter  the  orbit,  the  algorithm  switches  modes  and  then  tries  to 
maintain  equal  angular  spacing  during  orbit.  It  should  be  noted  that  this  is  just  one 
of  many  possible  solutions  to  cooperative  surveillance,  and  will  serve  as  a  baseline  for 
evaluating  performance  of  alternative  algorithms. 


Surveillance 

Orbit 


Figure  3.9:  Dubins  Path 
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3.6.2  Cooperative  Dubins  Paths.  The  Dubins  path  is  composed  of  no  more 
than  three  parts:  two  circular  arcs  and  a  straight  path.  In  order  to  make  the  paths 
of  all  UAVs  equal,  a  distance  is  needed  for  the  sum  of  all  three  parts.  From  Figure 
3.9,  the  path  uses  two  clockwise(CW)  circular  arcs  of  equal  radii(R)  at  each  end  of 
the  straight  section  (S).  The  path  length  (PI)  is 

Pl  =  S  +  A1R  +  A2R  (3.13) 

where  A\  and  A2  are  the  angles  swept  (in  radians)  at  each  end.  All  paths  for  each  UAV 
for  the  algorithm  will  be  based  on  this  path.  There  are  other  ways  to  construct  this 
type  of  path  by  altering  the  rotation  directions  at  each  end  and  also  having  differing 
radii  at  each  end.  Equal  radii  is  used  for  simplicity.  The  reason  for  CW  rotations 
has  to  do  with  collision  avoidance.  The  surveillance  orbit  must  be  counter  clockwise 
(CCW)  due  to  the  left  facing  side  camera.  As  multiple  UAVs  enter  the  surveillance 
orbit  CW,  they  are  forced  to  approach  from  the  outside  of  the  orbit,  minimizing  the 
chance  of  intersecting  paths.  No  other  method  of  collision  avoidance  was  done  as  a 
part  of  this  effort,  although  it  is  being  investigated  separately  at  AFIT. 

Creating  equal  path  lengths  is  done  by  altering  the  circular  arcs.  The  airframe’s 
maximum  effort  turn  determines  the  minimum  radius  that  can  be  used  in  the  path. 
Using  this  radius,  pathlengths  for  all  UAVs  are  calculated.  The  UAV  with  the  largest 
pathlength  becomes  the  baseline.  For  all  other  UAVs,  the  radius  is  incrementally  in¬ 
creased  until  all  paths  are  equal.  In  some  instances,  there  is  not  enough  angular  sweep 
to  increase  the  path  length  to  the  desired  length  given  the  positions  and  headings.  In 
this  case,  the  algorithm  will  add  a  complete  circle  to  one  end  and  adjust  the  radius 
until  the  desired  value  is  achieved. 

Figure  3.10  is  a  quick  illustration  of  the  end  product.  This  depicts  4  UAVs 
initially  orbiting  a  racetrack  on  the  right  (in  blue)  then  converging  on  a  target  orbit 
on  the  left.  The  algorithm  sends  each  UAV  a  set  of  GPS  waypoints  (red  X’s)  based 
on  a  path  equal  in  length  to  the  furthest  UAV.  Note  that  the  UAV  furthest  from  the 
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Figure  3.10:  Waypoints  based  on  Equal  Path  Lengths 

target  is  on  the  lower  right  and  has  the  smallest  radius.  Note  that  radii  increase  if 
the  UAV  is  closer  to  the  target.  Also  note  that  a  full  revolution  is  added  to  the  upper 
right  UAV  due  to  lack  of  available  angular  sweep.  The  final  headings  at  the  target  are 
determined  by  breaking  the  surveillance  orbit  into  equal  parts,  in  this  case  every  90 
degrees.  By  keeping  and  maintaining  the  surveillance  orbit,  any  angle  of  the  target  is 
kept  in  view  at  all  times. 

Simultaneous  arrival  simplifies  many  aspects  of  the  approach  to  the  surveillance 
orbit.  First,  it  places  the  UAVs  in  the  needed  position  as  they  arrive.  This  prevents 
the  need  for  large  adjustments  after  arrival  and  aids  in  collision  avoidance.  It  also 
provides  a  nice  simple  approach  to  path  generation  and  velocity  feedback.  If  the  paths 
are  all  equal,  and  all  UAVs  travel  at  the  same  ground  speed,  then  arrival  should  be  the 
same.  Usually  the  real  situation  alters  these  conditions,  so  velocity  feedback  is  used. 
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This  is  an  arbitrary  constraint,  and  simultaneous  arrival  is  not  the  only  approach, 
but  this  approach  is  simple  and  easy  to  implement.  Next,  the  details  of  creating  the 
paths  is  described. 

Procerus  Technologies  Kestrel  Autopilot  and  the  accompanying  VC  software 
use  GPS  to  exchange  position  telemetry,  thus  arises  the  need  to  translate  Cartesian 
coordinates  to  Latitude  (La)  and  Longitude  (Lo),  and  vice  versa.  Using  the  World 
Geodetic  System  1984  (WGS-84)  ellipsoid  [13]  and  a  small  angle  approximation,  a 
change  in  position  can  be  converted  from  La/Lo/h  into  Pn,Pe,Pd  using  [13] 


pN  =  (Rrn  +  h)ALa 
Pe  —  ( Rp  +  h)  cos  LaALo 
pD  =  Ah 


where 
Rm 
Rn  = 


a( 1  —  e2) 


(1  —  e2  sin2  La)3/2 
a 

(1  —  e2  sin2  La)3/2 


(3.14) 


a  is  the  semi-major  axis  of  the  Earth,  e  is  the  eccentricity  of  the  WGS-84  ellipsoid. 
This  equation  set  is  valid  only  for  distances  that  are  “small”  compared  to  the  size  of 
the  Earth.  Rearranging  the  above  and  solving  for  A  La  and  A  Lo,  the  conversion  from 
Cartesian  to  La/Lo/h  is 


La 

Lo 

h 


La0  + 
Lo0  + 


A  pN 

Rm  h 
A  pE 

(Rp  +  h)  cos  La0 


hQ  +  A  pD 


(3.15) 


This  equation  requires  that  the  La/Lo/h  of  the  Cartesian  origin  is  known  (La0j Lo0/h0). 

The  heading  angle  (CW  starting  from  due  north)  also  needs  a  conversion  to 
conform  with  the  standard  mathematical  conventions  (CCW  starting  from  the  x- 
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axis).  Let  the  *  denote  the  converted  heading  angle, 
conversions  to  and  from  are 


7T 

X  =  2  -  X 

_  7T 

X  =  g  “  X 


With  angles  in  radians,  the 


(3.16) 


With  the  position  in  Cartesian  coordinates  and  heading  angle  using  standard 
mathematical  conventions,  the  first  step  to  creating  the  path  is  finding  the  centers  of 
the  end  circles  (See  Fig  3.9)  given  the  initial  and  final  positions/headings.  Subscript  1 
denotes  quantities  in  the  initial  circle,  and  Subscript  2  denotes  quantities  in  the  final 
circle.  Since  the  paths  run  CW,  the  center  of  the  circle  will  always  lie  to  the  right  of 
the  UAV,  from  the  pilot’s  perspective.  Consequently  the  angle  from  the  UAV  to  the 
center  of  the  circle  will  be  y*  =  Xi  2  —  tt/2.  If  the  UAV  positions  are  (pei,2,Pni,2), 
the  positions  of  the  circle  centers  are 


(  PEcl,2  I  _  /m,2  +  -RC0s(Xii2  -7T/2) 

\Pnc1,2  j  \Pni,2  +  R  sin(xi,2  -  tt/2) 

With  the  circle  centers,  x.s  can  be  found.  The  heading  in  the  straight  section  is 
parallel  to  the  line  passing  through  the  circle  centers  since  the  radii  are  equal.  The 
straight  path  heading  is 


X*s  =  arctan 


/  P Nc2  ~  Pnc  1 

\  Pec'2  —  PecA 


From  this,  Al  and  A2  can  be  calculated: 


(3.18) 


Al  =  xl  ~  Xs 
A2  =  x’s  -  V 


(3.19) 


Al  and  A2  must  always  be  positive  to  generate  a  positive  length.  When  a  negative 
quantity  occurs,  27t  must  be  added  to  A1/A2. 
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The  position  of  the  two  tangency  points  is  needed  to  find  the  length  of  S.  Using 
the  circle  center  positions,  the  tangency  points  are  located  at 


I  PEtl,2 
\Pmi,  2 


PEc  1,2  +  Rcos(xs*i,2  +  tt/2) 
,PNci,2  +  Rsm(xsU  +  71 /2) 


Then  the  straight  section  distance  is 


(3.20) 


S  =  V(Pm2  ~  Pmi)2  +  (PEt2  ~  PEti)2  (3.21) 


Some  simplifying  geometry  occurs  when  S>2R.  In  this  case,  the  distance  between 
circle  centers  is  the  same  as  S  eliminating  the  need  to  find  the  tangency  points. 

Now  Equation  3.22  can  provide  the  total  path  length.  The  next  step  is  to  find 
the  appropriate  R  for  each  UAV  such  that  all  path  lengths  are  equal.  Finding  the 
necessary  R’s  is  an  iterative  process.  Once  the  path  lengths  are  calculated  using  the 
minimum  radius,  the  largest  path  length  ( Phase )  becomes  the  baseline  for  the  rest. 
Figure  3.11  illustrates  the  loop  used  to  find  the  correct  R  for  each  UAV.  As  the  top 
decision  block  implies,  the  path  lengths  are  not  exactly  the  same,  but  are  found  within 
a  defined  tolerance.  After  the  path  length  is  calculated,  the  algorithm  checks  to  see 
whether  that  value  has  exceeded  the  baseline.  If  so,  dr  changes  direction  and  is  halved. 
If  Pli  has  not  exceeded  the  baseline,  R  is  incremented  by  dr.  The  next  decision  block 
checks  whether  the  solution  is  converging.  Non-convergence  occurs  when  there  is  not 
enough  angular  sweep  in  the  given  path.  In  this  case,  R  is  reset  to  the  minimum 
radius  and  a  full  circle  is  added  to  the  initial  arc  in  the  path.  This  approach  is 
not  mathematically  elegant,  but  yielded  stable  results  under  all  circumstances.  The 
algorithm  repeats  until  all  UAVs  have  the  same  path  length. 

With  a  solution  for  all  desired  paths,  the  information  is  converted  back  into 
GPS  format  for  the  autopilot.  The  autopilot  needs  a  discrete  set  of  GPS  waypoints. 
This  algorithm  is  configured  to  provide  a  waypoint  for  every  90  degrees  of  sweep  in 
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the  arcs,  in  addition  to  a  waypoint  at  each  end  of  the  straight  section.  With  all  paths 
defined,  the  next  item  to  address  are  the  UAV  velocities. 


3.6.3  Handling  Wind  with  Feedback.  For  all  UAVs  to  arrive  into  the  surveil¬ 
lance  orbit  at  the  same  time  and  phased  equally  around  the  orbit,  the  ground  velocity 
(Vg)  must  also  be  the  same  for  all.  In  the  presence  of  wind,  maintaining  this  ground 
velocity  depends  on  the  heading  of  the  UAV.  To  find  the  necessary  airspeed  (W), 
vector  addition  states  that 

Vg  =Va  +  V\ y  (3.22) 


Using  the  above,  the  commanded  values  for  Vac  and  x°*  become 


Xc*  =  arctan 


Vg  sin  X*g  ~  Vw  sin  Xw 


Vc  = 

v  a 


COS 


Vg  cos  x*g  -  Vw  cos  Xw 
Vgcc°sx*g  ~  Vw  cos  Xw 

arctan 

V  Vg  cos  Xg-Vw  cos  Xw  )  ) 


(3.23) 


The  Kestrel  autopilot  provides  a  wind  estimate  in  the  telemetry  from  each  UAV. 
The  commanded  ground  velocity  will  come  from  cooperative  control  algorithm  (Eqns 
3.33  and  3.36).  Each  time  the  algorithm  issues  a  ground  speed  command,  it  will  be 
transformed  into  an  airspeed  command  using  the  above  equation  before  forwarding  to 
the  UAV  autopilot.  However,  the  feedback  is  actually  done  in  terms  of  ground  speed. 
During  flight  testing,  the  autopilot  commands  the  heading,  so  the  first  equation  in 
3.32  is  not  used. 


The  feedback  routine  that  maintains  ground  speed  while  approaching  the  surveil¬ 
lance  orbit  does  not  use  time  directly  to  accomplish  simultaneous  arrival.  By  using 
the  mean  path  length  ( Plr )  remaining  for  all  UAVs  from  the  orbit,  it  creates  a  floating 
reference  for  proportional  feedback.  The  commanded  ground  speed  for  the  ith  UAV 
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becomes 


nom 


kl{Plr  —  Plri) 


where 


(3.24) 


where  Plri  is  the  remaining  path  length  for  the  Ith  UAV,  and  k\  is  a  proportional 


gain.  Finding  the  remaining  path  length  uses  the  waypoints.  Figure  3.12  illustrates 
how  to  find  this  distance.  When  the  path  is  calculated,  an  array  is  created  that  stores 
the  remaining  distance  in  the  path  at  each  waypoint  (WP).  For  example,  d§  is  the 
distance  from  WP  5  to  the  surveillance  orbit,  and  dtowP5  is  the  straight  line  distance 
from  the  UAV  to  WP  5.  Then  Plr  =  dtowP5  +  d$.  A  whole  new  Dubins  path  is  not 
created  when  the  UAV  is  off  the  desired  path  or  is  not  at  the  needed  ground  speed, 
but  uses  the  initial  path  and  adjusts  the  heading  and  airspeed  accordingly. 

Once  the  UAVs  reach  the  surveillance  orbit,  we  need  to  find  the  correct  geometry 
of  this  orbit  to  maximize  the  sensor  orientation  and  FOV.  Please  refer  to  Figure 
3.13.  Assuming  a  coordinated  turn  and  constant  altitude,  the  lift  of  the  airplane 
provides  both  the  centripetal  acceleration  {Vg/r)  of  the  orbit  and  counters  gravity. 
Trigonometry  also  relates  the  bank  angle,  sensor  depression  angle  (y^),  and  altitude(h) 
to  the  orbit  radius  (r).  This  creates  two  equations,  with  two  unknowns  (r,  0) 


V2 

—  =  g  tan  0 


(3.25) 


r 


r 


tan(yd  +  0)  =  t 


If  you  assume  that  0  is  small,  then  after  a  little  algebra  a  quadratic  equation  with 
respect  to  r  appears.  Taking  the  solution  that  maximizes  r  (and  minimizes  0),  the 
solution  to  r  is 


r  = 


1 

2 


Yl)2 

9  ) 


(3.26) 
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This  result  will  be  the  orbit  radius  used  in  the  surveillance  orbit.  A  more  rigorous 
treatment  of  this  orbit  is  provided  in  [7]. 

Once  the  UAVs  enter  the  surveillance  orbit,  the  next  feedback  scheme  tries  to 
maintain  the  angular  spacing  of  the  orbit  with  respect  to  the  target  position.  Figure 
3.14  illustrates  the  algorithm.  One  of  the  vehicles  is  designated  as  the  reference  vehicle 
(UAV  1).  The  orbit  is  divided  into  equal  pieces  with  respect  to  UAV  1,  and  the  ith 
UAV  is  assigned  a  reference  angle  Arefi  (in  this  case  Aref  is  every  90  degrees  for  4 
UAVs).  The  commanded  ground  speed  for  the  ith  UAV  in  this  mode  is 

{Uiom  T  k'2^Arefi  Aj)  \Areji  Ai\  tol 

(3.27) 

Pnom  | Arefi  Ai\  U  tol 

Note  that  the  commanded  speed  is  centered  around  a  nominal  velocity,  and  that  there 
is  a  buffer  zone  around  a  tolerance.  For  the  case  of  UAV  3,  the  commanded  velocity 
would  be  greater  than  the  nominal  since  the  angle  difference  Areji  —  Aj  is  greater 
than  the  tolerance.  Conversely,  UAV  4  would  be  commanded  a  velocity  lower  than 
nominal. 

In  a  flight  testing  situation,  the  Kestrel  Autopilot  takes  care  of  certain  functions. 
When  given  a  set  of  waypoints,  the  autopilot  controls  the  commanded  heading  so  that 
the  vehicle  always  heads  toward  the  current  WP.  The  lower  level  control  loops  dictate 
the  control  surface  positions  and  attitude,  so  typically  there  is  no  need  to  provide  servo 
or  roll/pitch/yaw  commands.  The  commands  primarily  left  to  the  user  are  waypoints, 
loiter  points,  altitude,  and  airspeed.  As  the  two  feedback  equations  suggest,  airspeed 
control  is  a  central  part  of  the  algorithm. 
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Figure  3.11:  Equate  Pathlengths  Algorithm 
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Figure  3.12:  Calculating  Remaining  Path  Length 
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Figure  3.13:  Relating  Bank  Angle  to  Orbit  Radius 
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3.6.4  Software  Overview.  The  CC  algorithm  software  tries  to  maximize  the 
higher  level  capabilities  of  the  Kestrel  autopilot.  The  three  basic  capabilities  used 
are  waypoint  sets,  loiter  points,  and  airspeed  control.  As  stated  above,  waypoints 
provide  the  capability  to  follow  a  defined  path.  Loiter  commands  create  a  circular 
orbit  about  a  given  GPS  location  at  a  given  altitude.  Airspeed  control  gets  to  and 
holds  a  commanded  airspeed.  All  of  these  commands  (and  more)  are  contained  in  the 
PC  based  VC  software  that  comes  with  the  Kestrel  Autopilot  [1]. 

The  cooperative  control  (CC)  algorithm  will  be  a  program  that  runs  simultane¬ 
ously  with  VC.  Since  VC  communication  with  both  the  UAVs  and  external  software 
is  TCP/IP  based,  all  commands  are  network  packets  sent  to  the  UAVs  via  VC.  Also, 
all  telemetry  packets  from  the  UAVs  are  forwarded  from  VC  to  the  CC  program  using 
network  loopback.  The  CC  software  is  written  in  C++,  but  utilizes  custom  header 
hies  from  Procerus  that  enable  communication  with  VC,  and  also  special  header  hies 
from  Mathworks  that  enable  MATLAB  routines  to  be  utilized  in  C++.  Using  Visual 
Studio  2005  for  development,  the  software  is  a  Windows  based  and  event  driven  with 
a  graphical  user  interface  (GUI). 

The  MATLAB  capabilities  handled  the  complex  calculations  needed  to  create 
the  hnal  positions  and  headings  in  the  Dubins  paths,  to  create  the  waypoint  sets 
and  provide  airspeed  feedback  both  approaching  the  surveillance  orbit  and  in  the 
surveillance  orbit.  The  following  functions  were  written  as  a  part  of  this  effort  and 
are  used  in  the  software: 

•  mlfEnterOrbit  -  Using  the  target  location,  nominal  airspeed,  and  desired  alti¬ 
tude,  this  routine  calculates  the  orbit  radius  based  on  Equation  3.35  and  returns 
the  GPS  coordinates  and  headings  of  the  endpoints  of  the  Dubins  paths  that 
are  tangent  to  the  surveillance  orbit. 

•  mlfCreateOrbit  -  Using  the  current  UAV  locations  and  the  GPS  coordinates/- 
headings  entering  the  surveillance  orbit,  this  routine  returns  a  set  of  GPS  way- 


50 


points  with  Dubins  paths  of  equal  length.  The  algorithm  shown  in  Figure  3.12 
generates  the  WPs. 

•  mlfUpdateGrndSpd  -  Using  the  current  position  of  the  UAV  and  the  current 
waypoint  the  UAV  is  headed  to,  this  routine  returns  a  set  of  ground  speeds 
for  each  UAV  to  all  enter  the  surveillance  orbit  simultaneously.  Equation  3.33 
provides  the  ground  speeds. 

•  mlfAirspeed  -  Using  the  results  of  either  mlfUpdateGrndSpd  or  mlfMaintainOr- 
bit  and  the  current  wind  estimate,  returns  a  set  of  commanded  airspeeds  using 
Equation  3.32. 

•  mlfMaintainOrbit  -  Using  the  current  UAV  locations  and  target  location,  this 
routine  generates  a  set  of  commanded  ground  speeds  using  Equation  3.36. 

The  accompanying  CD  contains  the  C++  source  code  and  Appendix  B  contains  the 
supporting  MATLAB  routines.  Since  the  CC  software  is  event  driven,  four  events 
drive  this  program:  pushing  the  “Calc  Final  Points”  button,  pushing  the  “Upload 
to  VC”  button,  an  arrival  of  a  telemetry  packet,  and  a  timer  interrupt.  Figure  3.15 
illustrates  the  sequence  of  operations  for  each  event.  The  first  event  occurs  when  the 
“Calc  Final  Points”  button  is  pressed.  This  utilizes  the  mlfEnter Orbit  function  to 
create  endpoints  for  the  Dubins  paths  at  the  surveillance  orbit.  The  second  event, 
push  the  “Upload  to  VC”  button,  creates  WPs  based  on  equal  paths,  checks  them 
against  a  predefined  box  of  limits,  and  begins  the  process  of  uploading  all  WPs  with 
the  UploadCommands()  subroutine.  This  routine  controls  the  three  distinct  modes 
of  the  program  (See  Figure  3.16):  Generating  and  Uploading  WPs  (Mode  1),  arriving 
simultaneously  (Mode  2),  and  maintaining  orbit  (Mode  3).  The  sole  purpose  is  to 
govern  the  transition  from  one  mode  to  another.  The  reason  for  the  first  mode  is 
due  to  the  Kestrel  Autopilot.  It  will  not  start  executing  the  WP  set  until  the  whole 
set  is  uploaded.  WPs  are  only  uploaded  one  at  a  time,  and  the  next  WP  is  not  sent 
until  an  “Acknowledge”  (ACK)  packet  is  received  back  from  the  UAV  autopilot.  Once 
all  WPs  are  uploaded  to  all  UAVs,  the  program  transitions  to  mode  2.  Airspeed 
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Figure  3.15:  Software  Events 

commands  are  sent  every  second  to  aid  simultaneous  arrival.  Once  all  UAVs  arrive  on 
the  surveillance  orbit,  then  the  program  switches  to  mode  3  which  adjusts  airspeed 
every  second  based  on  angular  position  with  respect  to  the  target. 
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Figure  3.16:  Logic  Flow  of  UploadCommandsQ 


Figure  3.17:  CC  Algorithm  GUI 
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Referring  back  to  Figure  3.15,  the  third  event  (incoming  telemetry  packet)  is  key 
to  the  success  of  the  CC  algorithm.  An  incoming  telemetry  packet (TM)  causes  the 
program  to  update  the  state  variables.  Also,  when  the  program  sees  an  ACK  packet 
from  a  waypoint  command,  then  UpdateCommandsQ  can  send  the  next  waypoint  up. 
The  fourth  and  final  event  occurs  in  modes  2  and  3.  Once  mode  1  completes,  a  timer 
is  set  for  every  other  second.  Based  on  the  state  of  all  UAVs,  commands  are  sent  to 
adjust  the  airspeed.  The  only  way  the  program  can  transition  from  mode  3  back  to 
mode  1  is  to  push  the  “Upload  to  VC”  button.  If  the  user  so  desires,  he/she  can  enter 
a  new  target  and  the  UAVs  will  begin  the  algorithm  again.  A  sample  screen  shot  of 
the  GUI  is  shown  in  Figure  3.17. 

The  user  must  enter  certain  parameters  into  the  GUI.  In  the  “Target  Inputs” 
group,  the  program  needs  target  position,  surveillance  orbit  altitude,  desired  ground 
speed,  and  an  offset  angle  that  alters  the  position  of  first  UAV  position  in  the  orbit. 
An  offset  angle  of  zero  makes  the  first  UAV  position  due  East  of  the  target,  so  changing 
this  angle  rotates  this  position  counter-clockwise  about  the  target.  The  “Calc  Final 
Points”  button  uses  this  target  information  and  populates  the  final  lat/lon/heading 
for  each  UAV.  The  user  also  inputs  the  maximum  and  minimum  allowed  latitude/- 
longitude  for  waypoints  in  the  “Lat/Lon  Max/Min”  Group.  The  “Speed  Control” 
group  contains  the  parameters  for  mode  2  velocity  feedback,  and  the  “Orbit  control” 
group  in  the  GUI  contains  the  parameters  for  mode  3.  As  telemetry  is  received,  UAV 
information  is  automatically  updated  and  also  the  wind  information. 

The  CD  accompanying  this  thesis  contains  the  C-h-h  source  code  of  the  CC 
software. 

3. 7  Summary 

This  chapter  covered  the  approach  and  analysis  of  modelling  the  UAV  and 
creating  a  CC  algorithm  for  surveillance.  The  chapter  began  with  an  overview  of  both 
the  BATCAM  airframe,  the  Kestrel  Autopilot,  and  the  flight  testing  setup.  The  next 
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section  modified  a  first-order  model  from  Kingston  [9]  to  create  a  closed-loop  model 
of  the  BATCAM/autopilot  system.  Using  a  Dubins  based  path,  the  CC  algorithm 
provides  a  method  to  have  all  UAVs  arrive  at  the  surveillance  orbit  simultaneously, 
and  then  maintain  the  angular  spacing  during  the  orbit.  Software  written  in  C++ 
to  interface  to  the  VC  software  utilizes  the  analytic  power  of  MATLAB  to  create  a 
real-time  solution  to  implement  the  CC  algorithm.  The  next  chapter  will  summarize 
the  BATCAM  model,  in  addition  to  how  well  this  approach  worked  both  in  simulation 
and  in  flight  testing. 
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IV.  Results 


4-1  Introduction 

This  chapter  presents  the  results  of  modelling  and  the  performance  of  the  al¬ 
gorithm.  Using  actual  flight  test  data,  time  constants  for  the  first  order  responses 
are  extracted  and  woven  into  Equation  3.4.  Using  Simulink,  the  Cooperative  Control 
algorithm  is  connected  to  the  model  and  compared  to  flight  test  data  under  different 
wind  conditions.  Finally,  the  achieved  coverage  of  the  target  is  analyzed  for  both 
simulation  and  flight  test  results. 

4-2  Airframe  Model  and  Performance 

The  best  set  of  data  to  extract  model  parameters  occurred  on  the  last  flight 
of  testing  on  September  2,  2008.  The  time  of  day  was  near  sunset,  and  the  winds 
lessened  to  their  lowest  point  that  whole  day  (<  5  knots).  Flight  telemetry  contained 
less  disturbances  than  previous  flights  that  day,  and  was  considered  the  cleanest  pre¬ 
sentation  of  control  inputs  to  UAV  response. 

4-2.1  Closed-Loop  Property  Results.  The  model  contains  three  first-order 
time  constants:  roll,  pitch,  and  airspeed.  Using  the  time  lag  between  control  input 
and  response,  a  time  constant  can  be  extracted.  The  plots  of  control  input  and 
response  are  not  clean  step  responses.  Figures  4.1  and  4.2  show  a  small  section  of 
the  roll  and  pitch  commands,  respectively,  and  the  airframe  response.  Both  plots 
display  some  periodicity,  and  also  display  a  phase  difference  between  command  and 
response.  These  are  low  level  control  loops,  so  all  commands  originate  in  the  on-board 
autopilot.  This  is  important  because  without  ground  station  involvement,  the  delay 
between  when  the  command  is  issued  and  when  it  is  executed  is  “small”  compared  to 
the  period  shown  in  the  figures.  Assuming  negligible  delay,  from  basic  control  theory 
the  phase  lag  (Ph)  due  to  a  first-order  system  subject  to  a  periodic  forcing  function 
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Ph  =  —  arctan  a ;T 


is  [11] 


m  tan  Ph 

T  =  - 

co 


(4.1) 


where  co  is  the  frequency  of  the  input,  and  T  is  the  time  constant  of  the  first  order 
system.  The  plots  contain  both  Ph  and  co,  so  the  k  values  in  the  model  are  available 
by  finding  T.  co  is  calculated  finding  the  time  between  peaks  of  either  command  signal 
or  response  signal,  and  Ph  is  calculated  measuring  the  time  between  command  peak 
and  response  peak  using  the  period  of  co.  For  Roll  (Fig  4.1),  the  period  of  the 
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commanded  signal  was  2.3  seconds,  and  the  response  peak  lagged  the  command  peak 
by  1.2  seconds.  This  yielded  T=0.434  and  2.3  (from  Eqn  3.4).  For  pitch,  Figure 
4.2  yielded  a  command  period  of  10s  and  a  lag  time  between  command  and  response  of 
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Figure  4.1:  Typical  Roll  Response  to  Command 
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time  (s) 

Figure  4.2:  Typical  Pitch  Response  to  Command 

Is,  yielding  T=1.15  and  ko— 0.86.  Airspeed  commands  were  step  functions,  so  finding 
ky  was  more  straightforward(Figure  4.3).  Using  the  fact  that  reaching  95%  of  a  step 
function  takes  3T  for  a  first-order  system,  the  typical  rise  time  to  the  95%  level  took 
2.3  seconds  yielding  ky— 1.3.  This  was  true  for  both  large  (as  shown)  and  small  steps 
in  airspeed  command. 

This  approach  has  an  inherent  error  due  to  the  first  order  assumption.  Ex¬ 
amining  Figure  4.3,  the  response  appears  more  second  order  than  first  order.  The 
largest  consequence  is  that  the  model  will  not  capture  the  overshoot,  but  will  reach 
the  commanded  airspeed  similarly  to  the  actual  BATCAM/Kestrel  autopilot  system. 

The  model  also  contains  constants  (Section  3.4.2)  that  cannot  be  extracted 
from  flight  data,  but  are  either  chosen  by  the  user  or  are  determined  through  a  more 
subjective  tuning  process.  Maximums  and  minimums  on  commanded  roll  and  pitch 
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Figure  4.3:  Typical  Airspeed  Response  to  Command 

are  set  by  the  user  and  usually  based  on  common  practice.  Behavior  dealing  with 
tracking  the  intended  heading  or  altitude  are  based  on  user  preferences.  The  following 
values  were  chosen  based  on  creating  “reasonable”  behavior  in  the  model. 

The  roll  and  pitch  maximums  (Eqns  3.7  and  3.8)  are  set  to  ±30  degrees.  These 
were  the  limits  used  during  flight  test  and  were  intended  to  prevent  the  airframe 
from  reaching  an  uncontrollable  state.  One  very  practical  reason  for  preventing  an 
inversion  was  the  autopilot.  At  large  angles,  the  gyros  in  the  Kestrel  no  longer  yielded 
correct  attitude  measurements,  and  usually  caused  the  autopilot  to  act  erratically. 
When  angles  exceeded  45  degrees,  the  attitude  changes  very  quickly  and  verged  on 
instability.  30  degrees  seemed  a  reasonable  bound. 

An  upper  limit  to  airspeed  is  necessary  to  create  a  realistic  model.  If  you  exam¬ 
ine  Eqn  3.4,  without  an  upper  limit,  the  UAV  would  have  infinite  control  authority 
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with  respect  to  airspeed.  After  a  survey  of  flight  data,  Vramax=21.75  m/s  appeared  to 
be  the  maximum  value  with  full  throttle  without  being  in  a  dive  condition. 

The  next  values  ( kx,kh,a,XicPt )  were  determined  subjectively.  The  values  were 
based  on  desired  behavior.  For  kx  (see  Eqn  3.7),  the  desired  behavior  would  be 
to  command  a  bank  angle  that  would  not  saturate  until  there  was  a  difference  of 
about  30-40  degrees  between  commanded  heading  and  actual  heading.  For  this,  kx= 1 
accomplished  the  desire.  For  kh  (see  Eqn  3.8),  the  desired  behavior  would  be  not  to 
saturate  the  commanded  pitch  until  there  was  about  10  meters  difference  between 
commanded  altitude  and  actual  altitude.  This  led  to  a  value  of  k^= 0.05.  The  cross 
distance  parameters  were  also  based  on  the  10  meter  criteria  (see  Eqns  3.7,  2.3).  By 
letting  Xicpt  be  7t/4,  a  value  of  a=0.5  made  a  approach  7t/4  at  ?/s=10  meters.  10 
meters  seemed  reasonable  because  the  30  foot  deviations  allowed  some  margin  before 
maximum  control  effort  was  applied. 

This  first-order  approximation  of  the  closed-loop  BATCAM/autopilot  system 
was  chosen  for  a  number  of  reasons.  A  good  aerodynamic  model  does  not  exist  for 
this  airframe  and  the  Kestrel  autopilot  is  not  fully  documented.  The  BATCAM  has 
some  instabilities,  and  developing  a  controller  from  the  ground  up  that  both  mimicked 
the  autopilot  would  take  some  time.  With  two  major  missing  portions  of  the  closed- 
loop  model  (the  plant  and  controller),  it  seemed  reasonable  to  develop  a  model  based 
on  observed  data. 

Table  4.1  summarizes  the  model  parameters. 
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Table  4.1:  BATCAM  Model  Parameters 


Parameter 

Value 

Equation 

k# 

2.3 

3.4 

ke 

0.865 

3.4 

kv 

1.3 

3.4 

&max 

7t/6 

3.7 

@ max 

7t/6 

3.8 

V amax 

21.75  m/s 

- 

kx 

1 

3.7 

kh 

0.05 

3.8 

a 

0.5 

2.3 

Xicpt 

7t/4 

2.3 

h 

0.12 

3.24 

k2 

8 

3.27 
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4-3  Algorithm 

The  two  modes  of  the  algorithm  (approaching  the  orbit-Mode  2,  maintaining 
the  orbit-Mode  3)  have  two  different  performance  criteria.  The  first  is  measured  with 
cross  distance  and  arrival  time,  the  second  is  measured  with  cross  distance  and  angular 
position.  The  target  visibility  Section  (4.4)  will  analyze  target  visibility  based  on  the 
UAV’s  position  in  the  orbit.  Simulation  and  Flight  Test  results  are  analyzed  using 
this  criteria.  In  simulation,  four  wind  cases  are  tested:  no  wind,  10%  of  nominal 
airspeed,  25%  of  nominal,  and  50%  of  nominal.  While  orbiting,  the  algorithm  is  also 
tested  by  removing  a  UAV,  and  adding  a  UAV  to  the  formation.  The  criteria  for  this 
will  be  time  to  reconfigure. 

4-3.1  Simulation  Performance.  All  simulations  utilized  Mathwork’s  Simulink 
software  to  propagate  the  state  based  on  the  model  in  the  previous  section.  To  create 
simulations  for  each  mode,  two  Simulink  models  were  created.  The  control  algorithms 
in  the  first  include  velocity  feedback  to  aid  simultaneous  arrival  into  the  orbit.  The 
second  model  provides  airspeed  feedback  based  on  angular  position  in  the  orbit.  Fig¬ 
ure  4.4  displays  the  overall  model  as  it  appears  in  Simulink,  and  Figure  4.5  shows  the 
individual  BATCAM  model.  Associated  MATLAB  Code  is  in  Appendix  C. 

Blocks  labelled  BC  are  the  model  with  a  controller  for  roll,  pitch,  airspeed, 
desired  heading,  and  waypoint  location.  The  cooperative  control  controller  (CC  Con¬ 
troller)  takes  UAV  position,  yaw,  and  wind  information  then  feeds  the  individual 
BATCAM  models  waypoint  location,  commanded  heading,  commanded  airspeed,  and 
commanded  altitude.  Most  of  Matlab  function  blocks  in  the  simulation  are  nothing 
more  than  taking  the  data  vector  and  extracting  the  desired  elements. 

The  settings  for  the  Dryden  wind  model  varied  slightly  for  each  wind  case. 
Within  the  block  parameters  menu  inside  Simulink  the  wingspan  was  set  to  the  BAT¬ 
CAM  wingspan  of  21  inches  (adjusted  to  meters).  The  turbulence  probability  was  set 
to  “light”  for  the  first  three  wind  cases  (0,  10,  25%),  and  set  to  “moderate”  for  the 
50%  case.  All  other  settings  were  the  default  settings. 
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CC  Input  1x12 


Figure  4.4:  Overall  Simulink  Model 

The  first  simulation  takes  all  4  UAVs  and  tests  the  ability  to  arrive  at  the  same 
point  (Mode  2)  under  the  4  wind  conditions  (Wind  Heading  180  deg,  Wind  Speeds 
of  0,  1,17,  2.94,  5.88  m/s).  These  wind  speeds  are  10%,  25%,  and  50%  of  the  cruise 
speed  (11.75  m/s).  The  Dryden  wind  model  will  induce  randomness  into  the  wind 
speed  and  direction.  Table  4.2  displays  the  initial  conditions.  The  intended  4  paths 
are  straight  lines  starting  at  (pAr,p^)  =  (300,0)(250,0)(200,0)(0,50)  and  ending  at  the 
target  location.  All  UAVs  are  placed  slightly  off  their  intended  paths,  and  will  need  to 
drop  to  the  commanded  altitude  (see  Table  4.2).  The  simulation  ran  for  25  seconds  to 
provide  enough  time  to  reach  the  target.  The  scenario  is  set  up  to  test  the  algorithm’s 
ability  to  converge  and  stay  on  the  intended  path,  and  also  arrive  at  the  target  at  the 
same  time. 

Figures  4.6  and  4.7  present  the  position  of  all  4  UAVs  under  minimum  and 
maximum  winds  for  a  simulation  run  of  25  seconds.  Oscillations  along  the  intended 
path  increase,  but  all  UAVs  have  enough  control  authority  to  maintain  the  intended 
path.  Figures  4.8  and  4.9  plot  airspeed  vs  time  for  wind  speeds  of  0  and  5.88  m/s, 
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State  1x8 


Figure  4.5:  Individual  BATCAM  Simulink  Model 
Table  4.2:  Mode  2  Simulation  Initial  Conditions 


Item 

Location  ( Pn,Pe ) 

Heading 

Initial  Va 

Initial  Alt 

Target 

(200,200) 

n/a 

n/a 

n/a 

UAV1 

(290,-10) 

0 

11.54 

60 

UAV2 

(220,10) 

0 

11.54 

60 

UAV3 

(200,0) 

0 

11.54 

60 

UAV4 

(20,20) 

0 

11.54 

60 

Commanded  Values 

- 

Varies 

11.75 

50 

respectively.  A  quick  comparison  of  the  plots  illustrates  that  maintaining  the  path 
requires  much  more  control  effort  than  the  “no  wind”  case.  Maximum  required  air¬ 
speeds  differed  by  more  than  33%  with  respect  to  the  no  wind  case.  Table  4.3  presents 
the  performance  of  Mode  2  (Arrive  Simultaneously)  of  the  algorithm. 

Arrival  times  for  all  wind  cases  did  not  differ  by  more  than  2  seconds.  In  the 
worst  case  (wind=5.88  m/s),  the  time  difference  between  first  and  last  was  1.7  seconds. 
The  increasing  control  effort  displayed  in  the  airspeed  figures  resulted  in  sooner  arrival 
times.  Average  cross  distances  seemed  largely  dependent  on  initial  conditions,  but 
did  slightly  increase  with  increasing  wind  velocities.  The  minimum  distances  to  the 
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Wind  =  0  m/s 


300 


250 

~  200 
_£Z 

i _ 

O 

^  150 

100 

50 


Figure  4.6: 

target  dmin  increased  as  much  as  400%  between  the  0  wind  case  and  the  worst  wind 
case.  Depending  on  the  UAV  sensor’s  FOV,  this  miss  distance  could  affect  target 
visibility. 

Overall,  given  the  extracted  capabilities  of  the  BATCAM,  the  simulation  pre¬ 
dicts  reasonable  performance  of  the  Mode  2  algorithm.  Both  arrival  times  and  cross 
track  distances  appear  favorable  to  positioning  the  UAV  to  view  the  target.  The 
impact  of  UAV  attitude  on  target  visibility  will  be  addressed  in  the  Sensor  Footprint 
results.  It’s  important  to  note  that  this  model  is  reducing  a  complex  device  like  an 
autopilot  down  to  a  handful  of  equations.  Depending  on  the  actual  implementation 
of  the  Kestrel  control  loops  present  at  any  given  time,  behavior  could  vary  greatly. 
This  will  be  investigated  during  flight  tests. 


Simulation-Mode  2- Wind  0  m/s-UAV  position 
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Wind  =  5.88  m/s 


Figure  4.7:  Simulation-Mode  2-Wind  5.88  m/s-UAV  position 


Table  4.3:  Velocity  Feedback  Simulation  Results 


Item 

UAV  1 

UAV  2 

UAV  3 

UAV  4 

t arrive  Olll/ S 

21.60 

20.90 

21.03 

21.37 

t arrive  1 . 1 7m/ S 

21.34 

20.70 

20.87 

21.37 

t arrive  2.94m/s 

20.75 

20.08 

20.30 

21.05 

t arrive  5.88m/s 

19.45 

18.53 

18.72 

20.24 

mean  ys  W=0m/s 

10.57 

2.05 

4.54 

11.60 

mean  ys  W=1.17m/s 

10.31 

2.48 

4.59 

11.82 

mean  ys  W=2.94m/s 

9.81 

2.98 

4.55 

12.03 

mean  ys  W=5.88m/s 

9.30 

3.75 

4.82 

12.18 

dmin  W=0m/s 

0.84 

0.11 

0.11 

0.51 

dmin  W=1.17m/s 

0.17 

1.05 

0.89 

0.09 

dmin  W=2.94m/s 

0.75 

2.11 

1.96 

0.43 

dmin  W=5.88m/s 

1.46 

2.64 

4.35 

3.65 
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Airspeed  (m/s) 


Wind  0  m/s 


Figure  4.8:  Simulation-Mode  2-Wind  0  m/s-UAV  Airspeed  vs  Time 
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Wind  5.88  m/s 


Figure  4.9:  Simulation-Mode  2-Wind  5.88  m/s-UAV  Airspeed  vs  Time 
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Table  4.4:  Mode  3  Simulation  Initial  Conditions 


Item 

Location  (pN,pE) 

Heading 

Initial  Va 

Initial  Alt 

Target 

(0,0) 

n/a 

n/a 

n/a 

UAV1 

(-10,  65) 

I—1 

o 

o 

11.54 

50 

UAV2 

(70,0) 

281° 

11.54 

50 

UAV3 

(0,-45) 

191° 

11.54 

50 

UAV4 

(-90,10) 

101° 

11.54 

50 

Commanded  Values 

- 

Varies 

11.75 

50 

The  Mode  3  (maintaining  orbit)  simulation  will  begin  with  UAV  positions  shown 
in  Table  4.4  and  again  with  a  wind  heading  of  180°.  The  two  measures  of  performance 
will  be  cross  track  distance  and  angular  error.  The  target  will  be  located  at  the  origin, 
and  the  nominal  orbit  radius  of  63  m  gives  an  altitude  of  50m  (from  Eqn  3.35).  The 
UAV  at  the  6  o’clock  position  is  given  a  large  cross  track  error  outside  the  desired 
orbit,  and  the  UAV  at  the  9  o’clock  position  is  given  a  cross  track  error  inside  the 
orbit.  The  simulation  is  run  for  40  seconds  so  that  all  UAVs  complete  one  orbit 
around  the  target.  The  remaining  portion  of  this  section  presents  the  angular  errors 
introduced  from  the  addition  and  deletion  of  a  UAV.  The  Sensor  Footprint  Section, 
which  follows,  presents  the  target  visibility  metric  derived  from  the  data. 

The  impact  of  wind  starts  to  reveal  itself  during  orbiting.  Figures  4.10,  4.11, 
and  4.12  display  the  UAV  positions  over  the  40  second  simulation  under  0%,  25%,  and 
50%  of  nominal  airspeed.  The  0  wind  plot  shows  that  all  UAVs  are  able  to  converge  to 
the  desired  orbit  and  remain  there  with  little  error  in  crosstrack  ( ys ).  At  25%  winds, 
the  crosstrack  error  grows  but  all  UAVs  are  able  to  remain  on  orbit.  At  50%  winds, 
instabilities  start  to  occur.  UAV  4  is  never  able  to  recover  from  the  initial  large  initial 
crosstrack  error  (UAV  4  starts  at  the  6  o’clock  position).  On  the  upwind  side  of  the 
orbit  (right  side)  ys  decreases,  but  on  the  downwind  side  the  UAV  is  unable  to  stay 
on  orbit.  Instabilities  also  occur  with  UAV  3  (starting  at  9  o’clock  pos)  which  also 
began  with  a  large  crosstrack  error.  Table  4.5  presents  the  mean  angular  errors  (AE), 
mean  crosstrack  errors  (ys),  and  settling  time  (ts)  for  all  Mode  3  simulation  runs. 
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Settling  time  is  defined  as  the  time  required  to  converge  to,  and  remain  at  less  than 
±5  meters  ys  and  less  than  10°  AE.  From  Table  4.5,  none  of  the  UAVs  were  able  to 
settle  at  50%  winds.  Remember  that  the  desired  angular  position  is  based  on  UAV 
1,  so  there  is  no  angular  error  for  this  vehicle. 


Figure  4.10:  Simulation-Mode  3-Wind  0  m/s-UAV  position 

Figures  4.13,  4.14,  4.15,  and  4.16  illustrate  the  crosstrack  and  angular  errors  as 
functions  of  time.  Much  like  the  data  illustrated  in  the  positional  plots,  the  crosstrack 
plots  illustrate  the  effect  of  wind  on  the  ability  to  stay  on  orbit.  UAV  4  error  increases 
dramatically  on  the  downwind  side  nearing  40  meters  error.  During  that  same  section 
the  angular  error  was  negative,  forcing  a  command  to  increase  throttle  to  get  back 
into  position.  This  is  a  situation  where  the  algorithm  can  have  difficulties.  The 
downwind  side  of  the  orbit  increases  the  airspeed,  the  large  negative  angular  error 
will  also  force  a  further  increase  in  airspeed,  resulting  in  an  overshoot  that  may  be 
unrecoverable.  The  combination  of  high  winds  (50%  of  Vnom )  and  large  crosstrack 
errors  causes  instability  in  this  algorithm.  Winds  up  to  25%  of  Vnom  appear  to  be 
within  the  capabilities  of  the  BATCAM  and  the  algorithm.  If  the  UAV  starts  close 
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Wind  =  2.94  m/s 


Figure  4.11:  Simulation-Mode  3-Wind  2.94  m/s-UAV  position 


to  the  desired  position  (like  UAV  1),  then  the  algorithm  can  keep  position  within  the 
box  defined  by  ts  most  of  the  orbit  at  50%  winds.  Combining  the  results  of  Figures 
4.12,  and  4.16,  the  onset  of  instability  appears  to  occur  with  AE  greater  than  25° 
under  these  wind  conditions. 
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Wind  =  5.88  m/s 
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Figure  4.12:  Simulation-Mode  3-Wind  5.88  m/s-UAV  position 


Wind  =  0  m/s 


Figure  4.13:  Simulation-Mode  3-Wind  0  m/s-Crosstrack  Error  vs  Time 
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Wind  =  5.88  m/s 


Figure  4.14:  Simulation-Mode  3-Wind  5.88  m/s- Crosstrack  Error  vs  Time 


Wind  =  0  m/s 


time  (s) 


Figure  4.15:  Simulation-Mode  3-Wind  0  m/s- Angular  Error  vs  Time 
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Wind  =  5.88  m/s 


Figure  4.16:  Simulation-Mode  3-Wind  5.88  m/s-Angular  Error  vs  Time 


Table  4.5:  Mode  3  Simulation  Results 


Item 

UAV  1 

UAV  2 

UAV  3 

UAV  4 

AE  W=0m/s 

- 

4.5° 

7.6° 

OO 

o 

AE  W=1.17m/s 

- 

o 

oo 

7.9° 

OO 

bo 

o 

AE  W=2.94m/s 

- 

6.2° 

00 

o 

12.9° 

AE  W=5.88m/s 

- 

12.8° 

14.1° 

26.0° 

mean  ys(m)  W=0m/s 

0.47 

1.42 

1.78 

4.86 

mean  ys( m)  W=1.17m/s 

1.00 

1.74 

2.43 

5.48 

mean  ys( m)  W=2.94m/s 

1.36 

1.81 

3.09 

7.11 

mean  ys(m)  W=5.88m/s 

2.49 

2.61 

5.82 

17.30 

ts( s)  W=0m/s 

0 

2.85 

5.40 

15.36 

ts( s)  W=1.17m/s 

0 

2.52 

5.80 

21.00 

ts( s)  W=2.94m/s 

0 

8.1 

9.6 

- 

ts( s)  W=5.88m/s 

- 

- 

- 

- 
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Table  4.6:  Mode  3  Simulation  Initial  Conditions-Remove  UAV 


Item 

Location  (pN,pE) 

Heading 

Initial  Va 

Initial  Alt 

Target 

(0,0) 

n/a 

n/a 

n/a 

UAV1 

(-10,  65) 

I—1 

o 

o 

11.54 

50 

UAV2 

(70,0) 

281° 

11.54 

50 

UAV3 

(0,-53) 

191° 

11.54 

45 

Commanded  Values 

- 

Varies 

11.75 

50 

The  next  item  to  test  in  Mode  3  is  the  removal  of  a  UAV  from  the  formation. 
For  this  simulation,  the  scenario  will  reduce  a  4  UAV  orbit  to  3  UAVs.  The  initial 
conditions  are  given  in  Table  4.6,  the  wind  speeds  are  the  same  with  the  same  heading 
of  180°.  The  same  Mode  3  measures  of  performance  are  used:  crosstrack  distance  (ys), 
angular  error  (AE),  and  settling  time  (adjusting  time).  The  same  criteria  for  ts  will 
be  used  as  before  (ys  <  ±5m,  AE  <  ±10°).  The  UAVs  are  given  smaller  crosstrack 
errors,  are  placed  at  the  3  o’clock,  12  o’clock,  and  9  o’clock  positions  and  will  need 
to  re-adjust  to  120°  apart. 


Wind  =  0  m/s 


Figure  4.17:  Simulation- Mode  3/Remove-Wind  0  m/s-UAV  position 
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Wind  =  2.94  m/s 


Figure  4.18:  Simulation-Mode  3/Remove-Wind  2.94  m/s-UAV  position 


Figures  4.17,  4.18,  and  4.19  displayed  much  the  same  results  as  the  previous 
Mode  3  data:  the  algorithm  can  converge  at  wind  speeds  up  to  25%  of  Vnom.  The 
largest  difference  is  the  large  overshoot  of  UAV  3.  This  is  due  to  the  large  angular  error 
initially,  driving  a  large  increase  in  airspeed  to  correct  the  error.  Figure  4.22  provides 
a  clearer  picture  of  the  angular  error,  where  the  amount  of  necessary  correction  at 
t=0  is  slightly  over  50  degrees.  Also  like  the  previous  data  set,  instabilities  creep  in 
at  5.88  m/s  (50%  of  Vnorn)(Yig  4.19)  displayed  by  large  crosstrack  errors.  Again,  the 
downwind  side  of  the  orbit  exacerbates  overcorrections. 

The  time  needed  to  get  the  UAV  into  the  correct  position  is  shown  in  Table 
4.7.  The  zero  wind  case  sets  the  bench  mark  where  UAVs  2  and  3  are  in  position 
by  t=8.6  and  t=22.1,  respectively.  For  UAV  2,  the  time  to  get  into  position  actually 
decreases  with  increasing  wind.  The  geometry  of  the  problem  aids  this.  With  UAV 
2  starting  at  the  12  o’clock  position,  it  immediately  starts  turning  downwind,  so  the 
wind  actually  helps  the  UAV  get  into  the  desired  position.  From  Figure  4.23,  it 
actually  is  in  the  desired  position  in  2.9  seconds  and  remains  in  position  until  about 
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Wind  =  5.88  m/s 


Figure  4.19:  Simulation-Mode  3/Remove-Wind  5.88  m/s-UAV  position 

t=19s  when  UAV  1  (the  reference)  starts  down  the  downwind  side  of  the  orbit.  From 
an  angular  perspective,  UAV  3  does  not  perform  too  bad  (Fig  4.23).  It  reaches  the 
desired  position  at  t=25.5s  and  remains  within  10  degrees  of  desired  through  the  end 
of  the  simulation.  Unfortunately  at  this  time,  the  crosstrack  error  starts  to  exceed 
25m  (Figure  4.21).  Both  this  data  set  and  the  previous  display  that  the  BATCAM 
does  not  have  enough  control  authority  to  remain  in  position  when  winds  are  at  the 
5.88  m/s  using  the  Dryden  wind  model. 
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Wind  =  0  m/s 


Figure  4.20: 


Simulation-Mode  3/Remove-Wind  0  m/s-Crosstrack  Error  vs  Time 


Table  4.7:  Mode  3  Simulation  Results- Remove  UAV 


Item 

UAV  1 

UAV  2 

UAV  3 

mean  ?/s(m)  W=0m/s 

0.47 

1.57 

5.66 

mean  ys(m)  W=1.17m/s 

1.01 

1.75 

6.03 

mean  ys( m)  W=2.94m/s 

1.36 

1.70 

6.70 

mean  ys(m)  W=5.88m/s 

2.49 

2.60 

13.25 

ts( s)  W=0m/s 

- 

8.6 

22.1 

ts( s)  W=1.17m/s 

- 

7.9 

21.6 

ts( s)  W=2.94m/s 

- 

6.0 

25.5 

ts( s)  W=5.88m/s 

- 

- 

- 
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Wind  =  5.88  m/s 


Figure  4.21:  Simulation-Mode  3/Remove-Wind  5.88  m/s-Crosstrack  Error  vs  Time 


Wind  =  0  m/s 


Figure  4.22:  Simulation-Mode  3/Remove-Wind  0  m/s-Angular  Error  vs  Time 


79 


Figure  4.23:  Simulation-Mode  3/Remove- Wind  5.88  m/s-Angular  Error  vs  Time 
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Table  4.8:  Mode  3  Simulation  Initial  Conditions-Add  UAV 


Item 

Location  ( pN,pE ) 

Heading 

Initial  Va 

Initial  Alt 

Target 

(0,0) 

n/a 

n/a 

n/a 

UAV1 

(-10,  65) 

I—1 

o 

o 

11.54 

50 

UAV2 

(57,-32) 

281° 

11.54 

50 

UAV3 

(-51,-32) 

191° 

11.54 

45 

UAV4 

(-60,-20) 

135° 

11.54 

50 

Commanded  Values 

- 

Varies 

11.75 

50 

The  final  simulation  for  this  section  involves  inserting  1  UAV  into  a  3  UAV 
formation  and  observing  the  time  to  adjust  into  an  equally  spaced  4  UAV  orbit.  Initial 
conditions  are  shown  in  Table  4.8,  UAVs  are  placed  at  3,  11,  and  two  at  7  o’clock 
positions.  This  configuration  forces  two  of  the  UAVs  to  slow  down  to  approach  the 
9  and  12  o’clock  positions(with  respect  to  UAV1  at  3  o’clock),  and  the  final  UAV  to 
speed  up  to  reach  the  6  o’clock  position.  The  wind  conditions  will  be  the  same  as 
before,  and  small  errors  in  both  cross  track  and  heading  are  built  in.  Measures  of 
performance  will  again  be  mean  crosstrack  error  ys,  and  time  to  adjust  ts,  with  the 
same  settling  time  criteria  ( ys  <  ±5m,  AE  <  ±10°).  Angular  error  of  UAVs  2,  3,  and 
4  will  be  presented  as  a  function  of  time.  The  simulation  will  run  for  40  seconds. 

Figures  4.24,  4.25,  and  4.26  display  the  UAV  positions  for  the  simulation.  Typ¬ 
ical  for  past  simulations,  large  overshoots  of  the  orbit  do  not  begin  until  the  last  wind 
case  (|IU|=5.88  m/s).  UAV  4  is  commanded  to  increase  airspeed  to  get  to  the  desired 
position,  consequently  creating  an  overshoot  at  the  bottom  of  the  orbit  in  all  three 
Figures.  According  to  Table  4.9  this  UAV  does  reach  and  hold  the  desired  location  for 
the  first  two  wind  cases,  but  is  unable  for  the  last.  In  the  worst  case,  large  corrections 
in  heading  due  to  large  crosstrack  errors  again  lead  to  unstable  behavior.  For  the 
other  three  UAVs,  crosstrack  error  remains  fairly  tight  for  all  wind  cases,  with  mean 
values  under  6m. 

From  Table  4.9,  the  adjustment  times  remained  very  close  to  each  other  for 
UAVs  2  and  3  for  the  first  three  wind  cases,  but  UAV  4  displayed  increasing  ts  with 
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Figure  4.24:  Simulation-Mode  3/Add-Wind  0  m/s-UAV  position 

increasing  wind.  This  case  had  the  largest  initial  difference  in  angular  position  with 
UAV  3  approximately  66°  ahead  of  the  desired  position.  This  large  difference  did 
not  lead  to  increasing  adjustment  times  with  worsening  winds  (except  for  the  worst 
case).  It  appears  that  being  ahead  of  the  desired  position  is  more  favorable  than  being 
behind.  Increasing  airspeed  tends  to  lead  to  larger  overshoots  and  instability.  None 
of  the  UAVs  were  able  to  hold  their  position  within  the  designated  box  in  the  worst 
wind  case. 
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Wind  =  2.94  m/s 


Figure  4.25:  Simulation-Mode  3/Add-Wiud  2.94  m/s-UAV  position 


Table  4.9:  Mode  3  Simulation  Results-Add  UAV 


Item 

UAV  1 

UAV  2 

UAV  3 

UAV  4 

mean  ys(m)  W=0m/s 

0.61 

1.89 

1.27 

1.82 

mean  ys(m)  W=1.17m/s 

1.18 

2.24 

1.72 

2.41 

mean  ys{ m)  W=2.94m/s 

1.56 

2.17 

2.31 

3.91 

mean  ys( m)  W=5.88m/s 

2.70 

2.52 

5.75 

13.69 

ts( s)  W=0m/s 

- 

9.09 

14.16 

5.96 

ts( s)  W=1.17m/s 

- 

9.22 

14.38 

7.18 

ts( s)  W=2.94m/s 

- 

9.80 

14.29 

13.13 

ts( s)  W=5.88m/s 

- 

- 

- 

- 

83 


Wind  =  5.88  m/s 


Figure  4.26:  Simulation-Mode  3/Add-Wind  5.88  m/s-UAV  position 


Wind  =  0  m/s 


Figure  4.27:  Simulation- Mode  3/Add-Wind  0  m/s-Crosstrack  Error  vs  Time 
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Wind  =  5.88  m/s 


Figure  4.28:  Simulation-Mode  3/Add-Wind  5.88  m/s-Crosstrack  Error  vs  Time 


Wind  =  0  m/s 


Figure  4.29:  Simulation-Mode  3/Add-Wind  0  m/s-Angular  Error  vs  Time 
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Wind  =  5.88  m/s 


Figure  4.30:  Simulation-Mode  3/Add-Wind  5.88  m/s-Angular  Error  vs  Time 
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For  all  simulation  results,  UAVs  were  able  to  approach  the  orbit  for  all  wind 
cases,  but  were  unable  to  maintain  orbit  when  winds  were  50%  of  the  nominal  air¬ 
speed.  While  orbiting,  commanding  large  airspeeds  under  windy  conditions  led  to 
larger  crosstrack  errors  and  instability.  The  largest  crosstrack  errors  occurred  on  the 
downwind  side  of  the  orbit  and  led  to  instability  if  the  UAV  had  large  initial  crosstrack 
position.  The  algorithm  is  stable  and  can  maintain  orbits  when  winds  are  only  25% 
of  Vnom,  even  with  large  initial  errors.  Adjustment  times  for  a  63m  radius  orbit  when 
adding  and  deleting  UAVs  was  on  the  order  of  14-21  seconds  with  initial  angular 
errors  of  up  to  66  degrees. 

4-3.2  Flight  Test  Performance.  This  section  presents  the  position  accuracy 
of  the  algorithm  under  flight  test  conditions.  As  in  the  previous  section,  the  measures 
of  performance  are  crosstrack  error  (Mode  2  and  Mode  3),  arrival  time  (Mode  2),  and 
angular  error  (Mode  3).  The  algorithm  was  tested  two  times  in  flight  test.  The  first 
was  unsuccessful  because  airspeed  commands  from  the  algorithm  conflicted  with  the 
altitude  control  loop  in  the  autopilot.  Consequently,  algorithm  airspeed  commands 
were  ignored  due  to  the  tight  altitude  hold  parameters.  The  second  flight  test  success¬ 
fully  tested  two  BATCAMs  with  the  algorithm.  The  data  presented  in  this  section  is 
from  the  second  flight. 

Conditions  were  moderate  during  this  test.  Winds  were  out  of  the  East  South 
East  (Heading  of  300°)  with  speeds  of  3.0  to  3.6  m/s.  These  winds  were  30.6%  of  the 
steady  level  flight  speed  of  the  BATCAM.  Temperature  varied  from  40  to  50  degrees 
F,  skies  were  overcast,  and  visibility  was  approximately  1  mile  due  to  haze. 

Many  measures  were  taken  to  ensure  safety  and  positive  control  of  all  vehicles. 
UAVs  were  altitude  separated  by  100  feet  (30m)  as  a  collision  avoidance  measure.  Also 
the  CC  algorithm  software  was  outfitted  with  a  safety  measure  that  verihed  that  Lati¬ 
tude/Longitude  coordinates  for  each  WP  did  not  exceed  the  immediate  area.  Failsafe 
rally  points  were  established  so  that  the  BATCAM  would  loiter  a  predefined  point  if 
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communication  with  the  ground  station  was  lost.  The  UAVs  were  also  commanded 
to  begin  loitering  if  the  GPS  solution  was  lost. 


Figure  4.31:  Flight  Test  Location  at  Camp  Atterbury 


UAVs  were  flown  above  the  North  end  of  the  runway  keeping  all  UAVs  within 
vision  (See  Figure  4.31).  The  ground  station/trailer  was  parked  at  the  Southwest 
corner  of  the  pad  adjacent  to  the  airfield  operations  building.  The  distance  from 
ground  station  to  the  UAV  rarely  exceeded  1/4  mile.  The  GPS  limits  prevented 
overflight  of  the  road  to  the  north  and  east,  the  treeline  to  the  west,  and  the  line 
approximately  1/4  mile  south  of  the  trailer. 

Six  separate  tests  of  the  algorithm  were  conducted.  Two  tests  apiece  for  three 
different  airspeeds:  10  m/s,  11.75  m/s,  and  14  m/s.  The  first  flight  test  only  used  1 
nominal  airspeed  to  center  commands  around.  Analysis  from  this  first  flight  test  day 
implied  that  the  “slow  down”  commands  pushed  the  limits  of  the  BATCAM’s  ability 


to  maintain  lift.  The  three  separate  airspeeds  for  the  second  flight  test  day  tried  to 
quantify  if  algorithm  performance  was  affected  by  nominal  airspeed. 


Position:  2  BATCAMs,  t=325s  to  436s 


Figure  4.32:  Flight  Testing  -  Two  BATCAM  Test  1  -  Position 


Figures  4.32  and  4.33  present  the  position  plots  for  the  first  two  tests  with  two 
BATCAMs  (BCs).  Figures  4.34  and  4.40  illustrate  typical  orbit  position  plots  during 
flight  test.  Table  4.10  presents  the  average  measures  of  performance  for  Modes  2  and 
3.  The  six  tests  were  averaged  into  three  sets  of  results  based  upon  nominal  airspeed. 
As  a  comparison,  Table  4.11  displays  the  results  from  the  first  flight  test. 

As  stated  in  Chapter  III,  the  crosstrack  error  control  loop  (commanded  heading) 
resides  on  the  Kestrel  autopilot,  not  in  the  CC  algorithm.  Once  the  algorithm  provides 
the  autopilot  with  a  set  of  waypoints  for  all  UAVs,  the  only  commands  the  algorithm 
issues  are  airspeed  commands.  Internal  to  the  autopilot  is  the  method  to  keep  the  UAV 


Position:  2  BATCAMs,  t=934s  to  1054s 


Figure  4.33:  Flight  Testing  -  Two  BATCAM  Test  2  (Mode  2  and  3)  -  Position 


close  to  the  intended  path.  Unfortunately  a  detailed  description  of  the  commanded 
heading  method  was  not  provided  by  Procerus  for  the  Kestrel. 

Due  to  the  way  the  Kestrel  Autopilot  communicates  with  the  ground  station, 
it  took  approximately  7  to  9  seconds  for  the  VC  software  to  upload  the  waypoints  to 
the  BCs.  Consequently,  the  first  waypoint  often  appears  “misplaced”  since  the  CC 
algorithm  creates  a  path  off  the  most  recent  received  position.  During  the  delay,  the 
UAV  can  travel  as  far  as  100  meters.  This  is  due  to  the  wireless  TCP/IP  communi¬ 
cation  scheme.  Each  WP  is  a  packet,  so  if  the  algorithm  issues  7  WPs  for  each  UAV, 
the  communication  must  receive  an  “acknowledge”  back  from  the  autopilot  before 
sending  the  next  WP.  Also  the  autopilot  will  not  start  executing  the  flight  plan  until 
the  complete  plan  is  received. 
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The  Mode  2  results  for  the  two  BATCAM  set  differed  from  simulation  by  200 
to  300%.  In  the  simulation,  wind  results  for  25%  of  Vnom  as  a  comparison,  arrival 
times  for  all  UAVs  only  differed  by  1.6  seconds  (Table  4.3).  In  flight  test,  arrival  times 
differed  by  as  much  as  6  seconds.  Comparing  the  arrival  time  differences  for  the  three 
Thom’s,  the  ability  of  the  algorithm  to  get  UAVs  to  arrive  simultaneously  degrades  with 
increasing  Vnom.  This  trend  also  exists  in  the  crosstrack  error  results.  Just  like  linear 
control  theory  would  predict,  increased  velocity  results  in  larger  overshoots  when 
using  proportional  feedback.  Whether  the  desired  path  is  parallel  or  perpendicular 
to  the  wind  also  affects  crosstrack  error.  From  Figure  4.33,  the  long  straight  section 
of  BC  2’s  path  is  nearly  downwind,  leading  to  minimal  crosstrack  error,  but  BC  l’s 
straight  section  is  nearly  crosswind,  leading  to  larger  errors. 

Mode  3  results  showed  similar  results  to  Mode  2:  increasing  Vnom  increases 
positional  errors.  Angular  error  ranged  from  30°  to  40°,  and  crosstrack  error  ranged 
from  19  to  28  meters.  The  increased  operating  airspeed  tends  to  increase  the  overshoot 
distances.  The  shape  for  all  orbits  took  on  a  distinct  look. 

All  orbit  tracks  (see  Figures  4.40,  4.42,  4.44)  displayed  the  “egg  shaped”  orbit 
with  the  large  lobe  at  the  north  end  of  the  orbit.  Figure  4.34  contains  the  position  of 
both  BCs  for  1  orbit  and  a  plot  of  the  angular  error  as  a  function  of  BC  2’s  angular 
position  with  respect  to  the  target.  Both  BCs  are  orbiting  CCW.  BC  2  begins  the 
plot  with  a  headwind  at  the  7  o’clock  position,  with  BC  1  at  the  11  o’clock  position. 
The  angular  error  grows  for  the  first  half  of  the  orbit  because  BC  1  is  gaining  on  BC 
2.  As  BC  2  reaches  the  downwind  portion  of  the  orbit,  it  quickly  corrects  much  of  the 
angular  error  (AE)  with  the  superposition  of  increasing  airspeed  commands  to  correct 
the  AE  plus  the  wind  boost.  This  chain  of  events  appears  to  create  the  “egg.”  The 
target  visibility  section  (Section  4.4)  will  address  how  this  affects  target  visibility. 

In  general,  the  algorithm  produced  much  larger  crosstrack  errors,  angular  errors, 
and  arrival  time  differentials  than  simulation.  The  crosstrack  error  discrepancy  is 
likely  due  to  the  different  way  the  Kestrel  autopilot  corrects  crosstrack.  The  method  to 
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Table  4.10:  Flight  Test  2  Results-2  BATCAM  Test 


Vnom= 11-75  m/s 

Item 

BC  1 

BC  2 

mean  ys(m)  mode  2 

29.15 

20.03 

arrive  |  (s)  niod.6  2 

- 

5.05 

mean  AA(deg)  mode  3 

- 

32.33° 

mean  ys(m)  mode  3 

23.22 

18.83 

Vnom= 10  m/s 

Item 

BC  1 

BC  2 

mean  ys(m)  mode  2 

23.35 

17.71 

arrive  |  (s)  mod.6  2 

- 

1.94 

mean  A  A  (deg)  mode  3 

- 

30.7° 

mean  ys(m)  mode  3 

20.84 

18.82 

Vnam—14:  m/s 

Item 

BC  1 

BC  2 

mean  ys( m)  mode  2 

24.57 

22.21 

arrive |  (s)  mod.6  2 

- 

6.25 

mean  A  A  (deg)  mode  3 

- 

40.95° 

mean  ys(m)  mode  3 

28.40 

23.20 

correct  crosstrack  (Eqn  2.2)  in  the  simulation  probably  differs  from  Procerus’s  method 
programmed  in  the  Kestrel  autopilot.  Also  parameter  settings  can  affect  performance. 
Angular  errors  for  all  three  nominal  airspeeds  were  less  than  41°.  All  arrival  times 
were  within  6.25  seconds  of  each  other  with  the  least  difference  at  \4om=10rn/s. 


Table  4.11:  Flight  Test  1  Results-2  BATCAM  Test 


Item 

BC  2 

BC  3 

mean  ys( m)  mode  2 

10.45 

32.11 

|  arrive  \  (s)  mod.6  2 

- 

52.6 

mean  AA(deg)  mode  3 

- 

159.8 

mean  ys( m)  mode  3 

32.11 

23.64 
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Angular  Error  (deg)  North  (m) 


2  BATCAMs-Mode  3,Vo=12m/s,2nd  Test 


Angular  Error  vs  BC  2  Angular  position 


Figure  4.34:  Flight  Testing  -  Two  BATCAM  (Mode  3)-Typical  Orbit  -  Position 
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4-4  Target  Visibility 

The  final  section  of  this  chapter  presents  the  performance  of  the  algorithm  in 
terms  of  visibility  of  the  target.  In  addition  to  whether  the  target  is  in  the  FOV,  the 
target  will  also  only  be  visible  from  certain  directions.  For  this  analysis,  the  target 
will  always  be  located  at  the  origin  and  only  be  visible  for  headings  — 7t/2  to  tt/2  (from 
the  north).  The  measure  of  performance  will  be  the  time  these  conditions  are  met. 
For  simulation,  coverage  for  2,  3  and  4  UAV  formations  for  all  four  wind  conditions 
will  be  presented.  As  before,  errors  in  crosstrack  and  heading  are  introduced.  For 
flight  test  2  and  3  UAV  formations  were  used. 

4-4-1  Simulation  Results.  Two  approaches  will  be  presented  to  determine 
target  visibility.  The  first  will  be  a  more  restrictive  approach  that  determines  the  time 
the  target  is  visible  in  the  FOV  while  in  “the  box.”  Figure  4.35  defines  the  box  as  the 
angle  B  as  the  angle  from  the  vertical  where  2B=360/(Num  of  UAVs).  This  is  the 
angle  that  divides  the  surveillance  orbit  equally  based  on  UAV  quantity  and  decreases 
with  increasing  numbers  of  UAVs.  Three  box  angles  will  be  used:  B=90°,  B=60°,  and 
B=45°.  Ideally,  the  target  will  be  visible  the  whole  time  when  the  UAV  is  in  the  box. 
Two  metrics  will  be  measured,  time  in  the  box,  and  the  sum  of  all  times  for  each  UAV. 
When  compared  to  the  simulation  time,  this  will  provide  a  measure  of  effectiveness 
of  the  algorithm  compared  to  the  ideal  case.  The  second  approach  will  provide  the 
total  time  for  any  UAV  to  view  the  target  from  the  north  (B=90°),  regardless  of  the 
number  of  UAVs.  This  second  approach  (Overall  Total  time)  accounts  for  overlap 
time  between  any  two  UAVs.  From  a  human  factors  viewpoint,  the  first  approach 
keeps  the  ground  station  personnel  from  having  to  adjust  to  the  new  perspective  by 
keeping  1  UAV  as  the  primary  while  in  the  box.  The  second  approach  will  yield  higher 
in-view  times,  but  may  disorient  the  operator  with  camera  switching. 

Illustrating  the  motion  of  the  FOV  in  time  is  the  best  way  to  understand  the 
processes  that  affect  target  visibility.  Changes  in  velocity,  attitude,  altitude,  and 
heading  distort,  rotate,  and  move  the  FOV  leaving  the  in-view  time  metric  somewhat 
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Surveillance 

Orbit 


Figure  4.35:  The  Box  Angle  (B) 

lacking.  Unfortunately,  a  paper  document  is  somewhat  limited  in  its  capabilities  to 
present  this.  To  compensate  for  this,  X-Y  overhead  plots  that  show  a  sequence  of 
FOV  position  and  UAV  position/heading  will  be  used.  These  snapshot  sequences  will 
help  explain  the  phenomena  affecting  target  visibility. 

The  UAV  initial  conditions  for  each  simulation  are  the  same  as  used  in  Section 
4.3.  The  two  BATCAM  simulation  uses  the  same  initial  positions  as  UAV  1  and  UAV 
3  in  Table  4.4.  The  three  BATCAM  simulation  initial  conditions  are  shown  in  Table 

4.12.  The  four  BATCAM  case  also  uses  Table  4.4  for  beginning  values.  The  two  and 
four  vehicle  cases  will  follow  the  same  trajectories  as  shown  in  Figures  4.10,  4.11,  and 

4.12,  but  the  three  BATCAM  scenario  used  a  new  simulation.  The  induced  errors  in 
the  initial  conditions  give  a  comparison  of  how  the  crosstrack  and  angular  position 
control  loops  affect  target  visibility. 

Tables  4.13,  4.14,  and  4.15  present  the  in- view  times  and  percentages  for  the 
two  approaches.  Increasing  winds  dramatically  decreases  the  effectiveness  of  the  UAV 
formation  on  target  visibility.  For  all  UAV  formations,  the  target  was  visible  less  than 
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Table  4.12:  Sensor  Analysis-3  BATCAM  Initial  Conditions 


Item 

Location  ( pN,pE ) 

Heading 

Initial  Va 

Initial  Alt 

Target 

(0,0) 

n/a 

n/a 

n/a 

UAV1 

(-10,  65) 

I—1 

o 

o 

11.54 

50 

UAV2 

(57,-32) 

281° 

11.54 

50 

UAV3 

(-60,-20) 

135° 

11.54 

45 

Commanded  Values 

- 

Varies 

11.75 

50 

50%  of  the  time  in  the  “box.”  Increasing  the  number  of  UAVs  appears  to  make  the 
surveillance  more  robust  to  increasing  winds. 

The  overall  times  illustrate  some  interesting  trends.  Increasing  the  number  of 
UAVs  does  improve  the  total  target  visibility  times.  For  all  wind  cases,  the  highest 
percentages  lie  in  the  4  UAV  simulation.  The  increased  redundancy  aids  the  percent¬ 
ages,  but  not  by  large  amounts.  Looking  at  each  wind  case,  the  percentage  spread 
between  the  2  UAV  case  and  4  UAV  case  did  not  differ  more  than  10%.  It  appears 
that  the  extra  redundancy  aids  the  numbers  with  increasing  wind.  In  the  worst  wind 
case,  the  4  UAV  formation  had  10%  more  target  visibility  than  the  2  UAV  formation. 

Initial  conditions  plays  a  large  part  on  the  performance  of  the  in-view  time. 
Ideally,  the  zero  wind  cases  should  all  have  100%  visibility  for  all  40  seconds  of  sim¬ 
ulation.  In  the  two  BATCAM  0  Wind  case,  initial  conditions  caused  the  target  to 
not  be  in  the  FOV  for  2.34  seconds.  For  the  first  0.87  seconds,  UAV  1  was  not  in 
the  upper  half  of  the  circle,  and  UAV  2  was  proceeding  south.  As  UAV  1  entered 
the  upper  half,  from  snapshot  Fig  4.36  the  FOV  did  not  encompass  the  target  until 
t=1.74  seconds.  This  is  because  the  initial  conditions  that  put  the  UAV  outside  the 
orbit  and  roll  angle  of  zero  caused  the  control  loops  to  command  a  turn  rate  to  get 
back  on  course.  Consequently  the  turn  rate  increased  the  bank  angle  such  that  the 
FOV  went  too  far  East  to  see  the  target.  The  remaining  0.6  seconds  of  non-visibility 
comes  from  a  period  of  between  t=34.06  and  t=34.66  where  UAV  3  left  the  “box” 
before  UAV  1  entered  the  “box.” 
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Figure  4.36:  Sensor  Analysis-Two  BC  Orbit-BCl  Position  and  FOV 

Like  the  above  example,  crosstrack  error  has  a  dual  effect  on  target  visibility. 
When  the  UAV  is  not  on  the  orbit,  the  FOV  is  also  displaced.  The  command  to  correct 
the  crosstrack  error  makes  the  UAV  roll,  moving  the  FOV  again.  A  robust  control 
loop  for  controlling  crosstrack  error  has  a  negative  impact  on  target  visibility.  The 
right  balance  seems  to  keep  the  target  visible  for  small  errors,  and  gradual  corrections 
to  keep  positional  requirements. 

BC  3  for  the  3  UAV  formation  illustrates  the  affect  of  increasing  wind  on  the 
FOV  and  target  visibility.  Figures  4.37,  4.38,  and  4.39  illustrate  the  position  of  BC 
3  under  0,  2.94,  and  5.88  m/s  winds.  As  the  wind  produces  crosstrack  errors,  the 
frequency  of  correction  commands  increases.  These  manifest  as  roll  commands  that 
move  the  FOV  radially  outward  for  crosstrack  errors  outside  the  orbit,  and  radially 
inward  for  crosstrack  errors  inside  the  orbit.  For  a  radially  outward  movement,  the 
FOV  decreases  in  size,  and  conversely  a  radially  inward  movement  increases  its  size. 
Crosstrack  errors  outside  the  orbit  are  more  likely  to  move  the  FOV  so  that  the  target 
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is  not  visible.  As  the  figures  illustrate,  increasing  winds  tend  to  cause  outward  errors 
decreasing  target  visibility. 


Wind  =  Om/s 


Figure  4.37:  Sensor  Analysis-BC3  and  FOV,  Wind  0  m/s 


The  wind  direction  and  viewing  angle  were  chosen  to  accentuate  the  effect  of 
crab  angle  on  target  visibility.  The  largest  effect  appears  in  quadrant  II,  where  the 
UAV  first  experiences  the  downwind  side  of  the  orbit.  One  way  to  alleviate  the  effects 
of  crab  angle  is  to  change  the  shape  of  the  orbit  from  circular  to  elliptical  much 
like  Rysdyk  [16]  and  Farrell  [7].  Rysdyk  varied  the  radius  as  a  function  of  heading, 
altering  altitude.  Farrell  maintained  constant  altitude.  Both  approaches  came  up 
with  ellipse-like  solutions  and  would  improve  target  visibility,  but  do  not  address  an 
approach  to  handle  crosstrack  errors.  The  circular  approach  was  chosen  due  to  the 
Kestrel  autopilot.  The  command  set  in  the  autopilot  had  a  loiter  command,  and  to 
create  an  ellipse  would  require  construction  with  many  waypoints.  The  delays  already 
experienced  with  the  Dubins  path  would  be  compounded  with  10+  more  waypoints 
to  upload.  It  seems  feasible  to  merge  the  ideas,  but  time  did  not  allow  during  this 
research. 


Table  4.13:  Simulation  Results-Target  Visibility  Time  with  2  BATCAMs 


Wind  Condition 

BC  1(%) 

BC  3(%) 

Box  Total(%) 

Overall  Total(%) 

0  m/s 

21.39(53.5%) 

16.93(42.3%) 

38.32(95.8%) 

37.66(94.1%) 

1.17  m/s 

18.27(45.7%) 

16.17(40.4%) 

34.44(86.1%) 

33.75(84.4%) 

2.94  m/s 

17.66(44.2%) 

12.70(31.8%) 

30.36(75.9%) 

29.82(74.6%) 

5.88  m/s 

12.23(30.6%) 

5.12(12.8%) 

17.35(43.4%) 

17.35(43.4%) 

all  results  are  in  seconds,  percentages  (%)  are  with  respect  to  total  simulation  time  40s 

Table  4.14:  Simulation  Results-Target  Visibility  Time  with  3  BATCAMs 


Wind  Condition 

BC  1 

BC  2 

BC  3 

Box  Total (%) 

Overall  Total(%) 

0  m/s 

13.72 

11.24 

11.23 

36.19(90.5%) 

38.25(95.6%) 

1.17  m/s 

11.83 

10.71 

11.03 

33.57(83.9%) 

36.82(92.1%) 

2.94  m/s 

11.21 

9.22 

8.84 

29.27(73.1%) 

34.63(86.6%) 

5.88  m/s 

7.67 

3.41 

2.48 

13.56(33.9%) 

18.82(47.1%) 

all  results  are  in  seconds,  percentages  (%)  are  with  respect  to  total  simulation  time  40s 

Table  4.15:  Simulation  Results-Target  Visibility  Time  with  4  BATCAMs 


Wind  Cond. 

BC  1 

BC  2 

BC  3 

BC  4 

Box  Tot. 

Overall  Tot.(%) 

0  m/s 

9.54 

10.76 

8.47 

8.36 

37.14(92.9%) 

38.63(96.6%) 

1.17  m/s 

7.71 

10.42 

8.01 

8.00 

34.15(85.4%) 

37.58(94.0%) 

2.94  m/s 

7.09 

9.11 

7.12 

3.94 

27.26(68.2%) 

36.23(90.6%) 

5.88  m/s 

4.52 

9.18 

1.39 

2.67 

17.76(44.4%) 

26.01(65.0%) 

all  results  are  in  seconds,  percentages  (%)  are  with  respect  to  total  simulation  time  40s 
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Wind  =  2.94m/s 


Figure  4.38:  Sensor  Analysis-BC3  and  FOV,  Wind  2.94  m/s 

Overall,  roll  angle  lias  a  large  effect  on  target  visibility  with  a  camera  fixed  to 
the  fuselage.  Even  in  no  wind  cases,  initial  conditions  caused  the  target  to  leave  the 
FOV  due  to  changes  in  roll.  The  largest  impact  came  from  crosstrack  errors,  where 
the  control  loops  altered  roll  angle  to  correct  positional  error.  As  wind  increased,  crab 
angle  started  to  play  a  role  in  decreasing  visibility  in  the  portions  of  the  orbit  where 
a  cross  wind  was  prevalent. 
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Wind  =  5.88m/s 


Figure  4.39:  Sensor  Analysis-BC3  and  FOV,  Wind  5.88  m/s 
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Table  4.16:  Flight  Test  1  Results-Target  Visibility  Time 


Test 

BC  1 

BC  2 

BC  3 

Box  Total 

Overall  Total(%) 

Tot  Orbit  Time 

2  BC  Test  1 

- 

15.57 

10.98 

26.55(25%) 

24.10(22.7%) 

106.2 

2  BC  Test  2 

- 

13.03 

16.42 

29.45(25.6%) 

28.00(24.3%) 

115.2 

3  BC  Test  1 

10.01 

11.03 

13.65 

29.57(25.5%) 

35.64(30.8%) 

115.8 

all  resul 

s  are  in  seconds,  percentages  (%)  are  with  respect  to  total  orbit  time 

4-4-2  Flight  Test  Results.  This  section  takes  flight  test  data  and  determines 
the  amount  of  time  the  target  was  visible.  Like  the  simulation,  the  target  will  only  be 
visible  from  the  north,  so  when  the  target  to  UAV  heading  is  between  90°  and  270° 
the  observed  time  is  not  counted.  The  same  measures  of  performance  will  be  used: 
in- view  time  in  the  “box” ,  total  in-view  time  in  the  box,  and  overall  time.  The  “box” 
will  be  defined  as  B=90°  for  2  BC  tests  and  B=60°  for  3  BC  tests.  Overall  time  is 
for  any  BC  where  B=90°. 

Like  Section  4.3.2,  most  of  the  results  come  from  the  second  test  flight,  but  as  a 
comparison,  Table  4.16  displays  results  from  the  first  flight  test.  As  stated  previously, 
the  airspeed  feedback  failed  to  work  on  the  first  flight  test.  Another  big  difference  was 
the  proportional  gain  was  increased  from  2  to  8  to  match  simulation  parameters.  Keep 
in  mind  that  the  CC  algorithm  does  not  provide  heading  commands,  only  airspeed 
commands  to  the  autopilot.  The  Kestrel  manages  the  crosstrack  error  correction 
(commanded  heading) . 

Table  4.17  presents  the  average  visibility  times  and  percentages  for  the  three 
Vnom  cases,  Figures  4.40  through  4.45  present  the  orbit  tracks  and  FOV  snapshot 
plots.  All  of  the  orbit  tracks  displayed  the  “egg”  with  the  large  lobe  primarily  on  the 
north  end  of  the  orbit.  Figures  4.40,  4.42,  and  4.40  distinctly  show  a  growing  large 
lobe  with  increasing  Vnom.  With  the  “box”  placed  on  the  north  end  of  the  orbit,  this 
large  lobe  caused  decreased  target  visibility.  The  14m/s  plots  did  display  a  rotation 
of  the  large  lobe  in  the  direction  of  the  wind,  which  increased  target  visibility  for  that 
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run  since  it  brought  the  trajectory  closer  to  the  nominal  orbit.  In  general,  when  the 
lobe  was  oriented  north,  poor  visibility  times  resulted. 

Just  like  the  simulations,  the  roll  commands  necessary  to  correct  crosstrack 
errors  wreak  havoc  on  the  target  visibility  times  (See  FOV  snapshot  Figures  4.41,  4.43, 
and  4.45).  Large  crosstrack  errors  on  the  large  side  of  the  “egg”  create  commands 
that  bring  the  FOV  toward  the  UAV  position.  Consequently,  the  target  quickly  leaves 
the  FOV  until  the  roll  angle  decreases.  Also  the  climbs  to  correct  altitude  also  affect 
the  FOV.  Whenever  the  FOV  moves  in  the  direction  of  travel  and  appears  ahead  of 
the  UAV  position,  the  UAV  is  pitching  up.  Pitching  up  seems  most  prevalent  in  the 
transition  from  crosswind  to  headwind  (and  vice  versa)  conditions.  Sometimes  this 
helps  target  visibility,  sometimes  it  hurts  visibility. 

In  its  attempt  to  keep  angular  position  while  orbiting,  the  algorithm  negatively 
impacted  target  visibility.  Comparing  the  first  flight  test  data  (Table  4.16)  to  the 
second  flight  test  data  (Table  4.17),  very  poor  angular  positioning  still  yielded  slightly 
better  target  visibility.  One  possible  reason  is  that  the  proportional  gain  for  the 
second  flight  test  may  be  too  high.  Crosstrack  overshoots  due  to  increasing  airspeeds 
due  to  angular  error  caused  rolling  of  the  airframe,  impacting  visibility.  As  is,  the 
proportional  gain  used  in  Equation  3.27  that  worked  well  in  simulation  did  not  work 
well  in  flight  testing.  Unfortunately,  the  simulation  does  not  closely  represent  the 
behavior  of  the  Kestrel  in  crosstrack  correction,  so  either  the  Kestrel  needed  further 
tuning  to  minimize  crosstrack  error  or  more  testing  is  needed  to  find  the  gain  that 
better  balances  positional  feedback  with  target  visibility. 

Despite  the  non-functioning  airspeed  feedback,  the  first  flight  test  did  display 
some  results  that  correlated  with  simulations  (Table  4.16).  Since  the  first  test  flew 
both  2  BCs  and  3  BCs,  the  data  displayed  that  target  visibility  increased  with  increas¬ 
ing  number  of  UAVs.  2  BC  formations  had  target  visibility  times  around  23%  and 
3  BC  formations  increased  the  percentage  to  30%.  Due  to  the  large  angular  errors, 
“box”  totals  saw  no  net  benefit. 
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Overall  the  multi-UAV  system  had  poor  results.  For  the  second  flight  test,  box 
totals  never  exceeded  25%  of  the  orbit  time,  and  overall  totals  never  exceeded  23% 
(Table  4.17).  Without  correct  positional  placement  and  attitude,  the  FOV  cannot  see 
the  target.  Crosstrack  errors  have  the  largest  impact  on  target  visibility.  Not  only 
do  they  move  the  FOV  radially,  but  force  the  autopilot  to  roll  the  fuselage  to  correct 
for  the  error.  With  a  body  fixed  camera,  roll  commands  to  correct  crosstrack  errors 
decreases  in-view  times.  As  seen  in  the  first  flight  test,  increasing  numbers  of  UAVs 
does  provide  a  certain  amount  of  robustness  and  increases  target  visibility. 
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Table  4.17:  Flight  Test  2  Results-Target  Visibility  Time 


Test 

BC  1 

BC  2 

Box  Total 

Overall  Total(%) 

Tot  Orbit  Time 

Test  1  14iom=H-75m/s 

14.42 

20.54 

34.96(21.9%) 

32.44(22.7%) 

159.5 

Test  2  Vnom=11.75rri/s 

10.30 

10.81 

21.11(21.8%) 

20.54(21.2%) 

96.7 

Test  3  Vn0m = 1  Orn/s 

8.44 

11.03 

19.47(18.5%) 

19.61(18.6%) 

105.2 

Test  4  14om=l  Om/s 

12.81 

7.10 

19.90(20.8%) 

18.38(19.2%) 

95.6 

Test  5  Kom=14m/s 

7.25 

10.61 

17.86(15.2%) 

16.37(13.9%) 

117.4 

Test  6  14om=14m/s 

12.38 

19.50 

31.88(25.7%) 

27.71(22.3%) 

123.85 

all  results  are  in  seconds,  percentages  (%)  are  with  respect  to  total  orbit  time 

2  BATCAMs-Mode  3,Vo=12m/s,2nd  Test 


Figure  4.40:  Flight  Testing  -  Orbit  Traces  for  14om=11.75m/s  (Mode  3) 


4-4-2. 1  Flight  Test  2. 
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Vnom=1 1 ,75m/s,  2nd  Test-BC  2  FOV 


Figure  4.41:  Flight  Testing  -  FOV  Snapshot  for  Vn0m= 11.75m/s 


2  BATCAMs-Mode  3,Vo=10m/s,2nd  Test 


Figure  4.42:  Flight  Testing  -  Orbit  Traces  for  for  Vnom= 1  Orn / s 
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Vnom=10m/s,  2nd  Test-BC  2  FOV 


Figure  4.43:  Flight  Testing  -  FOV  Snapshot  for  14om=10m/s 


2  BATCAMs-Mode  3,Vo=14m/s,2nd  Test 


Figure  4.44:  Flight  Testing  -  Orbit  Traces  for  14om=14m/s 
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Vnom=14m/s,  2nd  Test-BC  2  FOV 


Figure  4.45:  Flight  Testing  -  FOV  Snapshot  for  Vlom=14m/s 
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4-5  Summary 

This  chapter  presented  the  results  of  the  first  order  model  of  the  BATCAM, 
chosen  parameters  for  the  CC  algorithm,  performance  of  the  CC  algorithm  in  sim¬ 
ulation,  and  also  flight  test  performance.  Performance  was  quantified  from  both  a 
positional  and  target  visibility  perspective. 

Using  basic  control  theory,  flight  test  data  yielded  first  order  constants  for  air¬ 
speed,  roll,  and  pitch.  Other  parameters  that  closed  the  loop  and  mimicked  the 
Kestrel  autopilot  behavior  were  derived  from  actual  limits  set  during  flight  testing, 
and  reasonable  desired  behavior. 

The  first  round  of  simulation  analyzed  the  ability  of  the  BATCAM  model  to 
correct  and  maintain  desired  position  under  wind  conditions  of  0%,  10%,  25%  and 
50%  of  nominal  airspeed.  Position  accuracy  was  based  on  arrival  time  and  mean 
crosstrack  error  for  Mode  2  (simultaneous  orbit  approach),  and  based  on  settling 
time,  mean  angular  error,  and  mean  crosstrack  error  for  Mode  3  (maintaining  orbit). 
For  Mode  2  simulation,  UAV  arrival  times  differed  by  less  than  0.7  seconds  for  the  no 
wind  case,  and  less  than  2  seconds  for  the  50%  wind  case.  Crosstrack  error  steadily 
increased  as  wind  increased,  but  was  largely  dependent  on  initial  conditions.  For  Mode 
3  simulation,  the  CC  algorithm  was  able  to  maintain  orbit  for  the  first  three  wind 
cases,  but  instabilities  appeared  for  the  50%Vnom  case.  Again,  mean  crosstrack  error 
and  mean  angular  error  showed  a  strong  dependence  on  initial  conditions.  Mean 
angular  error  for  all  UAVs  rose  from  6.9°  for  0  wind  to  17.6°  for  50%  wind,  mean 
crosstrack  error  rose  from  2.11  m  for  0  wind  to  7.05  m  for  50%  wind.  Settling  time  for 
the  4  UAV  formation  was  15.36  seconds  for  o  wind,  and  never  settled  into  the  proper 
positions  for  the  worst  wind  case  due  to  large  crosstrack  errors.  These  instabilities 
begin  with  the  combination  of  moderate  crosstrack  error  on  the  downwind  side  of  the 
orbit. 

Performance  of  the  CC  algorithm  was  flight  tested  under  three  nominal  air¬ 
speeds:  10  m/s,  11.75  m/s,  and  14  m/s.  Performance  in  both  Mode  2  and  Mode  3 
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showed  increasing  error  with  increasing  Vnom  where  average  arrival  time  differentials 
rose  from  1.94  to  6.25  s,  average  angular  error  rose  from  30.7°  to  40.9°,  average  Mode 
2  crosstrack  errors  rose  from  17  to  22  m,  and  average  Mode  3  crosstrack  errors  rose 
from  18  to  23  m.  Surveillance  orbits  took  on  an  “egg”  shaped  appearance  due  to 
cycle  of  increasing  angular  error  followed  by  an  angular  correction  when  wind  aided 
the  increased  commanded  airspeed. 

Target  visibility  simulations  analyzed  2,  3,  and  4  UAV  formations  under  all  4 
wind  conditions  and  quantified  time  when  the  target  was  in  the  FOV.  Two  approaches 
were  taken:  one  measuring  the  time  the  UAV  was  in  “the  box” ,  the  other  measuring 
the  time  any  UAV  could  see  the  target  from  the  north  side  of  the  orbit.  Overall 
target  in-view  times  ranged  from  96.6%  for  the  4  BC  zero  wind  case  to  43.4%  for 
the  2  BC  50%  wind  case.  Percentages  not  at  100%  for  the  0  wind  case  was  due  to 
initial  conditions.  Increasing  the  numbers  of  UAVs  improved  in-view  times  as  wind 
increased.  For  50%  winds,  4  UAV  formations  kept  the  target  visible  65%  of  the 
time  verse  only  43%  for  the  2  UAV  formation.  Crosstrack  error  had  a  dual  effect  on 
target  visibility.  Large  crosstrack  errors  outside  the  orbit  both  translated  the  FOV, 
decreased  its  size,  and  caused  radial  motion  away  from  the  target.  Roll  commands  to 
correct  the  crosstrack  error  pointed  the  FOV  more  directly  downward  shrinking  size 
and  moving  it  radially.  So  as  crosstrack  error  increased,  target  visibility  decreased. 

Large  crosstrack  errors  due  to  the  “egg  shaped”  orbits  during  flight  testing 
resulted  in  poor  target  visibility.  Much  like  simulations,  the  necessary  roll  to  correct 
crosstrack  errors  moves  the  FOV  away  from  the  target.  Since  the  large  lobe  of  the 
“egg”  resided  primarily  on  the  north  side  of  the  orbit,  target  visibility  times  ranged 
from  14  to  23%  of  the  total  orbit  time.  The  first  flight  test  did  display  increasing 
target  visibility  with  increased  number  of  UAVs,  much  like  simulations. 
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V.  Conclusion 


This  research  accomplished  three  major  tasks:  developed  a  closed  loop  model  for  the 
BATCAM,  developed  a  cooperative  control  algorithm  to  control  multiple  BATCAMs 
to  accomplish  surveillance,  and  analyzed  the  performance  of  the  algorithm  in  both 
simulation  and  flight  test  environments. 

5. 1  Conclusions 

The  algorithm  met  its  overall  goal  of  commanding  multiple  UAVs  in  real-time 
to  approach  and  maintain  an  orbit  about  a  fixed  target.  During  flight  testing  UAV 
arrival  times  at  the  surveillance  orbit  differed  by  an  average  value  of  4.4  seconds  and 
maintained  angular  position  while  orbiting  within  an  average  value  of  34°.  Wind 
conditions  averaged  30%  of  Vnom  for  these  results. 

Simulation  data  showed  mixed  results  in  predicting  flight  test  phenomena.  Sim¬ 
ulation  did  predict  the  increasing  orbit  overshoots  with  increasing  wind,  and  also 
predicted  decreased  target  in-view  time  with  these  overshoots.  Simulation  did  not 
predict  the  extent  of  crosstrack  error  that  existed  in  flight  test  data. 

Wind  conditions  had  a  large  impact  on  the  positional  performance  of  the  algo¬ 
rithm.  In  simulation,  both  mode  2  and  mode  3  showed  increasing  errors  in  position 
and  target  visibility.  In  mode  2  mean  crosstrack  errors  displayed  only  small  increases 
(4.4%)  from  0  wind  to  50%  wind,  but  mode  3  mean  crosstrack  errors  jumped  331% 
from  2.13m  to  7.05m.  For  the  worst  wind  case,  mode  3  never  converged  to  maintain 
position  in  orbit,  and  was  unable  to  re-adjust  the  UAVs  in  the  add  or  remove  UAV 
simulations.  The  algorithm  was  able  to  maintain  position  and  adjust  at  25%  winds. 

The  simulation  wind  cases  pointed  out  that  the  useable  wind  speed  range  for  the 
algorithm  lies  between  zero  and  at  least  25%  of  Vnom.  Somewhere  between  25%  and 
50%  of  Vnom  instabilities  begin  especially  in  mode  3  (maintain  orbit).  If  the  crosstrack 
errors  exceed  a  certain  amount,  then  the  algorithm  overshoots  the  intended  path  on 
the  downwind  side  of  the  surveillance  orbit  and  is  unable  to  maintain  the  desired 
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orbit.  Flight  test  data  also  displayed  these  overshoots  when  wind  averaged  37%  of 

V  nom- 

Both  simulation  and  flight  test  data  show  that  increasing  numbers  of  UAVs 
can  improve  target  in- view  percentages  with  increasing  winds.  Overall  percentages  in 
the  worst  wind  case  between  the  2  BC  simulation  and  the  4  BC  simulation  differed 
by  21.6%.  Also  increasing  winds  had  less  of  an  effect  with  more  UAVs.  When  winds 
increased  from  0  to  50%  of  Vnom,  the  four  UAV  formation  only  lost  31%  visibility  time 
where  the  2  UAV  formation  lost  51%  of  target  visibility  time.  The  “box”  approach 
showed  mixed  results  where  regardless  of  the  number  of  UAVs,  in-view  percentages 
remained  similar  for  all  wind  cases. 

Crosstrack  error  had  the  largest  impact  on  target  visibility  due  to  its  dual  effect 
on  the  FOV.  The  error  itself  translates  the  FOV  away  from  the  desired  point  and 
the  roll  command  to  the  error  induces  further  FOV  movement.  For  the  purposes  of 
keeping  the  target  in- view,  the  crosstrack  feedback  loop  must  be  designed  to  balance 
the  need  to  correct  crosstrack  error  and  roll  commands  that  cause  target  visibility 
loss. 

Roll  has  the  largest  affect  on  the  FOV  for  a  body  fixed  side  camera.  Since  roll 
not  only  translates  the  FOV  but  also  alters  the  size  of  the  FOV,  maintaining  roll  and 
minimizing  changes  in  roll  becomes  very  important.  Due  to  the  size  and  mass  of  the 
BATCAM,  this  can  be  challenging.  As  both  simulation  and  flight  test  shows,  roll 
caused  the  target  to  leave  the  FOV,  especially  with  increasing  winds. 

5.2  Recommendations  for  Future 

Expanding  on  this  work  could  go  in  many  directions.  One  direction  refines  this 
work,  the  next  creates  a  new  branch  of  research  with  the  same  basic  elements.  The 
specialized  nature  of  Cooperative  Control  means  many  branches  of  research  can  stem 
from  the  same  idea  by  just  altering  the  scenario,  operating  environment,  or  even  the 
assumptions. 
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This  research  left  some  unfinished  work.  A  full  6  DOF  model  was  started  and  lies 
in  Appendix  A,  but  work  remains  in  closing  the  loop  and  finding  a  stable  controller 
that  mimics  the  Kestrel.  The  approach  took  DeLuca’s  [5]  aerodynamic  work  and 
linearized  it  to  generate  the  aerodynamic  forces  and  moments.  Look  up  tables  is 
another  approach  for  these  forces  and  moments.  A  controller  remains  undesigned. 

Due  to  the  limits  of  flight  testing  at  Camp  Atterbury,  the  “magic  combination” 
of  control  loop  settings  was  never  found  that  maximized  target  in- view  time  using  the 
built  in  loiter  command.  Since  the  graduates  of  the  March  2009  group  were  the  first 
to  fly  the  BATCAM,  more  than  half  the  flight  testing  entailed  creating  a  stable  set  of 
gains  for  the  autopilot.  Work  remains  to  hone  these  gains  for  surveillance  purposes. 

Mode  2  is  based  on  the  arbitrary  constraint  of  simultaneous  arrival.  This  ap¬ 
proach  simplifies  certain  aspects  of  the  problem,  but  is  not  the  only  way  to  approaching 
the  surveillance  orbit.  One  possible  modification  of  the  algorithm  is  to  change  Mode  2 
such  that  the  UAVs  arrive  asynchronously.  Both  collision  avoidance  and  readjustment 
become  major  components  of  the  problem,  but  may  be  a  more  flexible  approach. 

An  optimal  solution  to  this  problem  will  provide  a  gauge  as  to  the  quality  of  this 
heuristic  approach.  This  comparison  can  quantify  the  upper  limit  to  this  concept  in 
terms  of  performance.  As  shown  in  Section  2.3.2  [15],  finding  the  right  simplifications 
to  the  optimal  problem  can  possibly  yield  an  approach  worthy  of  real-time  use. 

This  research  illustrated  that  there  is  a  trade  space  to  explore  in  minimizing 
crosstrack  error  and  maximizing  target  in-view  time.  The  current  algorithm  only 
strives  for  positional  accuracy,  and  does  not  account  for  its  affect  on  FOV  movement. 
Finding  this  balance  between  positional  accuracy  and  keeping  the  FOV  on  the  target 
given  the  constraints  of  the  hardware  (BATCAM/Kestrel)  will  take  some  work  both 
in  simulation  and  flight  testing. 

This  vein  of  research  is  highly  specialized  and  requires  pulling  together  knowl¬ 
edge  from  many  disciplines.  Hands  on  work  in  Cooperative  Control  can  be  quite 
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challenging,  yet  rewarding.  No  textbook  or  simulation  can  provide  the  experience  of 
integrating  software  with  hardware. 
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Appendix  A.  Unfinished  BAT  CAM  Analysis 


A.l  Mass  Properties 

Mass  properties  include  the  mass,  the  center  of  gravity  location,  and  the  three 
moments  of  inertia.  The  first  two  properties  are  measured  using  direct  measurements: 
one  from  a  scale,  and  the  second  from  a  simple  balance  test.  The  moments  of  inertia 
are  derived  from  angular  velocity  measurements. 

Figure  A.l  illustrates  the  equipment  used  to  measure  angular  velocity.  This  ro¬ 
tational  table  made  by  PASCO  comes  with  DataStudio  software  and  a  rotary  motion 
sensor  that  can  measure  the  angular  position,  velocity,  and  acceleration.  Measure¬ 
ments  are  taken  at  0.1  second  intervals.  The  table  is  accelerated  by  weights  attached 
to  string  that  rest  on  the  outer  pulley  and  wind  around  the  main  rotational  axis. 
Three  different  weights  are  used  to  accelerate  the  table:  500  grams,  200  grams,  and 
50  grams. 


Figure  A.l:  Moment  of  Inertia  Measurement  Equipment 

The  procedure  is  as  follows.  Once  the  BATCAM  is  mounted  on  the  table  in  the 
desired  orientation,  the  software  is  set  to  start  recording  angular  velocity.  The  weight 
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is  attached  to  the  twine  and  the  table  is  allowed  to  freely  rotate  for  a  few  seconds. 
The  software  graphs  the  data  and  displays  the  angular  velocity  as  a  function  of  time. 
Since  air  drag  begins  to  affect  measurements  at  high  speeds,  only  the  data  in  the  linear 
portion  of  the  graph  is  used.  The  software  calculates  a  linear  fit  on  the  selected  data 
and  the  slope  (angular  acceleration)  is  recorded.  The  BATCAM  is  removed  and  the 
steps  are  repeated  to  find  the  inertia  of  the  table  itself.  Three  runs  are  conducted  for 
each  weight  for  a  total  of  9  runs  per  axis.  The  final  angular  acceleration  (O50o, 200,50) 
is  the  mean  of  all  three  runs  for  each  weight.  The  tension(T)  in  the  twine  is 

-<500,200,50  —  \mw  —  mf){9 - 2 - /  (A.l) 


where  mu,  is  the  mass  of  the  weight,  rn f  is  an  equivalent  mass  to  account  for  friction, 
g  is  the  acceleration  of  gravity,  and  d  is  the  diameter  of  the  spindle  at  the  base  of  the 
table.  The  moment  of  inertia  (I)  becomes 


-^500,200,50 


<^500,200,50 

2^500,200,50 


(A.2) 


The  final  moment  of  inertia  of  the  BATCAM  for  each  axis  averages  the  individual 
MOIs  from  each  weight  less  the  MOI  of  the  table. 


hoo  +  L200  +  h 


'50 


Lxx,yy,zz 


‘■table 


(A.3) 


Realize  that  Itabie  is  calculated  in  the  same  manner  as  other  MOIs,  just  without  the 
BATCAM  on  the  table.  The  results  for  all  needed  dimensions  and  mass  properties 
are  shown  in  the  table  below. 


A.2  Aircraft  Modelling -Deluca 

Anthony  DeLuca  [5]  completed  a  thorough  aerodynamic  analysis  of  the  BAT¬ 
CAM  airframe.  I11  addition  to  the  Lift  and  Drag  data  for  various  airspeeds  and 
angles  of  attack,  the  research  also  included  are  the  control  surface  moments  for  dif- 
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Table  A.l:  Dimensions  and  Mass  Properties 


Item 

Value 

Wingspan  (b) 

21  in 

Mean  Geometric  Chord 

5.94  in 

Wing  Reference  Area 

103.7  in2 

Mass 

0.425  kg 

I XX 

0.0028  kg  m2 

lyy 

0.0088  kg  m2 

Izz 

0.0091  kg  m2 

Table  A. 2:  Summary  of  Stall  Angles 


10  mph  a staii 

20  mph  a staii 

30  mph  astaii 

50  mph  astaii 

Flex  Wing 

8.7° 

12.7° 

14.8° 

7.5° 

Rigid  Wing 

12.8° 

12.7° 

8.5° 

- 

ferent  elevon  deflections,  and  thrust  coefficients  for  the  propeller /motor.  Since  the 
BATCAM  wing  is  flexible,  much  of  the  analysis  focuses  on  the  differences  between 
this  wing  and  a  rigid  wing. 

Comparisons  between  rigid  and  flexible  wings  revealed  that  flexible  wings  delay 
stall  conditions  for  3  out  of  the  4  airspeeds  (20,  30,  50  mph).  When  aerodynamic 
forces  deflect  the  trailing  edge  upwards,  the  relative  angle  of  attack  (AoA)  decreases. 
The  exception  to  this  was  at  10  mph,  where  stall  occurred  at  8.7°  verse  12.8°  AoA 
for  a  rigid  wing.  Apparently  a  laminar  separation  bubble  manifests  itself  as  slight 
undulations  in  the  lift  line,  degrading  aerodynamic  efficiency.  Table  A. 2  summarizes 
the  AoA’s  where  stall  conditions  occurred. 

Figure  A. 2  illustrates  the  lift  and  drag  plots  for  the  BATCAM  at  the  above 
airspeeds  at  various  angles  of  attack(a).  The  50  mph  data  was  incomplete  due  to 
wing  damage  at  the  last  data  point.  Note  the  increased  drag  at  10  mph  due  to  the 
same  degraded  aerodynamic  efficiencies  as  stall. 

The  BATCAM’s  V-tail  introduces  some  coupling  of  control  moments  compared 
to  the  traditional  aileron,  elevator,  and  rudder.  Three  [5]  cases  were  examined:  de- 


117 


q 

c 

o 

u 


-15 


- Flex  10  mph 

—  Flex  20  mph 

- Flex  30  mph 

- Flex  50  mph 


20  25 


-0.5- 


Angle  of  Attack  (deg) 


0.7 


5 

c 

« 

o 

3; 

8 

o 

0> 

a 


- Flex  10  mph 

- Flex  20  mph 

- Rex  30  mph 

—  Flex  50  mph 


-15 


-10 


-5  0  5  10 

Angle  of  Attack  (deg) 


15 


20 


25 


Figure  A. 2:  BATCAM  Flexible  Wing  CL  vs  a  and  vs  a  [5] 
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Table  A. 3:  Average  Slopes  of  Elevon  Runs 


30  mph 

Single 

Tandem 

Opposed 

dCm/d5e 

0.0087 

0.0133 

0.0020 

dQ/dSe 

0.0012 

0.0005 

0.0021 

dCn/d5e 

-0.0043 

-0.0015 

-0.0072 

flection  of  a  single  elevon,  deflection  of  both  elevons  in  the  same  direction  (tandem), 
and  deflection  of  both  elevons  in  opposing  directions.  Ideally,  opposing  deflections 
create  rolling  and  yawing  moments  with  negligible  pitching  for  heading  and  direction 
changes.  Tandem  deflections  create  pitching  with  minimal  rolling  and  yawing  for 
attitude  and  pitch  control. 

For  the  single  deflection,  all  three  moment  directions  behaved  nearly  linearly  up 
to  5°.  In  both  tandem  and  opposing  cases,  linear  behavior  held  up  to  8  °.  Table  A. 3 
presents  the  average  slope  per  degree  of  deflection  (<5e)  for  the  non  dimensionalized 
moment  coefficients  of  roll(Q),  pitch  (Cm),  and  yaw  ( Cn ). 


Figure  A. 3:  BATCAM  Thrust  Coefficient  vs  Advance  Ratio  [5] 
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The  thrust  coefficients  (CV)  were  very  consistent,  linear,  and  well  behaved. 
Figure  A. 3  presents  Ct  at  20  mph.  The  power  coefficients  are  plotted  as  a  function 
of  the  advance  ratio  Ja,  which  is  defined  as: 

Ja=—1  (A. 4) 

7]d 

where  Va  is  the  airspeed,  y  is  the  prop  speed  in  revolutions  per  second,  and  d  is  the 
diameter  of  the  prop.  At  a  =  0  the  BATCAM  would  overcome  a  20  mph  head¬ 
wind  at  approximately  1/2  throttle  (8000  RPM).  Axial  force  measurement  showed  no 
dependence  on  angle  of  attack. 

DeLuca  presented  a  very  solid  aerodynamic  analysis  for  use  in  model  build¬ 
ing.  The  above  will  be  used  to  create  linear  relationships  of  aerodynamic  coefficients 
centered  around  level  flight  conditions.  From  these  relationships,  an  estimate  for 
aerodynamic  forces  (Lift,  Drag,  Thrust,  etc)  will  be  incorporated  into  the  model. 

A. 3  Linearized  Aerodynamic  Coefficients 

This  section  pulls  data  from  DeLuca  [5]  presented  in  Chapter  II  and  creates 
linear  equations  for  aerodynamic  coefficients  around  level  steady  flight.  These  non 
dimensional  coefficients  provide  the  means  to  find  the  forces  and  moments  on  the 
airframe.  The  major  aerodynamic  forces  and  moments  on  an  airframe  are  :  Lift(L), 
Drag(D),  Side  Force  (Y),  Moments  due  to  control  surface  deflections(L,  M,N),  and 
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Thrust (Th).  These  quantities  are  defined  [18]  [5]  as: 


L  =  qSCL 
D  =  qSCD 
Y  =  qSCY 

L  =  qSbCi  (A.  5) 

M  =  qScCM 
N  =  qSbCN 
Th  =  CTpr)2d4 


where  q  is  the  free  stream  dynamic  pressure(=pV)j2/2),  S  is  the  wing  reference  area,  b 
is  the  Wingspan,  c  is  the  mean  geometric  chord,  p  is  the  air  density,  r)  is  the  propeller 
speed  in  revolutions  per  second,  and  d  is  the  propeller  diameter. 


The  coefficients  will  be  linearized  with  their  respective  independent  variable 
(either  AoA,  Side  slip  angle,  Advance  Ratio,  or  elevon  deflection).  For  the  BATCAM, 
level  steady  flight  occurs  at 


Vao  =  11.75  m/s 
Throttle  =  50% 


(A.6) 


Using  this  data,  the  nearest  appropriate  data  from  DeLuca  is  chosen.  Since  11.5  is 
approximately  25  mph,  either  the  20mph  or  the  30  mph  can  be  used.  For  this  research, 
all  linearized  coefficient  equations  will  use  0  degrees  AoA  and  20  mph  curves  from 
Deluca’s  data  (Figure  A. 2).  Drawing  a  straight  line  from  the  points  at  0  deg  AoA  to 
5  deg  AoA,  the  linearized  Cl  and  Cp  as  functions  of  a  (in  radians)  are 


(a)  =  5.730a  +  0.7 
CD(a)  =  1.203a  +  0.075 


(A.7) 


From  Deluca  [5],  the  side  force  coefficient  (Cy)  remained  constant  in  the  AoA  range 
of  0  to  5  degrees.  Also  for  side  slip  angles  (/?)  up  through  12  degrees,  Cy  increased 
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linearly.  The  relationship  between  Cy  and  /3( in  radians)  was 


CY  =  -0.573 f3 


(A.8) 


For  the  thrust  coefficient,  from  Figure  A. 3  the  linearized  equation  as  a  function  of  Ja 
becomes 

Ct(Jcl)  =  -0.1  Ja  +  0.135  (A.9) 


Next  is  to  gain  equations  for  Cl,  Cm,  and  Cm  as  a  function  of  elevon  deflection. 
Using  the  derivatives  shown  in  Table  A. 3  the  coefficient  equations  become 


r  sc, 

Cl  =  W.s‘ 

n  —  r 
Cm  ~ 

r  9Cn  a 

Cn  ~  aT5' 


(A. 10) 


For  this  analysis,  the  elevons  will  only  be  used  in  some  linear  combination  of  tandem 
and  opposing  configurations.  Realize  that  the  partials  for  the  tandem  set  and  the 
opposing  set  must  be  kept  together  in  order  to  describe  the  coupling  effect  of  the 
V-tail.  In  theory  [5],  a  tandem  command  should  only  create  a  pitching  moment. 
DeLuca’s  data  showed  a  light  rolling  and  a  yaw.  For  the  sake  of  model  building 
dCL/d8e  and  dCM/dde  will  be  zero  to  eliminate  any  asymmetries  that  may  have 
existed  in  DeLuca’s  UAV.  Also  in  theory,  opposing  elevons  should  only  create  roll 
and  yaw.  Once  again  for  opposing  elevon  commands,  dCM/dSe  will  be  zero  in  the 
model.  So  for  tandem  elevon  deflections 


CL  =  05, 
CM  =  0.01335, 


(A. 11) 


and  for  opposing  deflections 


CL  =  0.002 5e 

CM  =  05e  (A.12) 

CN  =  0.00725e 

Typically  any  small  eccentricities  in  control  surfaces  or  the  airframe  are  trimmed  out 
to  get  the  UAV  to  fly  straight. 

This  is  a  limited  set  of  the  many  actual  aerodynamic  forces  and  moments,  but 
does  capture  the  major  components.  Forces  and  moments  will  use  these  coefficient 
linearizations  and  apply  equation  3.5  for  calculation. 

A. 4  6  DOF  Model 

The  6  DOF  flat  earth  equations  [18]  are 

Fb 

vB  =  —0BvB  +  BBg0  H - 

m 

ujB  =  —J  lVtBJuB  +  J  1Tb 

(A. 13) 

Phi  =  E($)ub 

PNED  =  BbVB 

where  vB,  u>B  are  the  velocity  and  angular  velocities  in  the  body  frame,  $  is  the 
attitude  vector  (roll(</>),  pitch(6)),  yaw(-0)),  and  Pned  is  the  position  vector  in  the 
NED  frame.  FB  and  TB  are  the  body  forces  and  body  torques,  which  are  usually 
aerodynamic  ( FB  =(T-D  Y  -L)r,  TB  ={L  M  N)T).  Bb  is  the  transpose  of  C£  from 
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Eqn  3.11.  The  other  associated  matrices  are  (c=cos,  s=sin,  t=tan): 


Ixx 

0 

-i 

J  = 

0 

lyy 

0 

-Ixz  0 

4 

0 

-R 

Q 

= 

R 

0 

-p 

rQ 

P 

0 

1  tOscf) 

tdccj) 

0 

C(j) 

—  S(ft 

n 

S(j) 

C(f) 

VJ 

cO 

c6 

(A. 14) 


where  uB  =  (P  Q  R)t  and  the  moments  of  inertia  are  from  the  mass  properties  section. 
Ixz  was  not  measured,  but  is  usually  small  compared  to  the  other  three  moments  and 
may  be  set  to  zero  or  a  percentage  of  Ixx. 

To  properly  calculate  aerodynamic  forces  and  moments  use 


vR  =  vB  -  Bt 


WN 

WE 

WD 


(A. 15) 


relative  velocity  vr. 

This  is  the  open  loop  model.  The  controller  and  closed  loop  model  are  unfin¬ 
ished.  The  MATLAB  implementations  of  the  aerodynamic  forces  and  the  flat  earth 
6  DOF  equations  are  in  the  next  section. 
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A. 5  MATLAB  files 


Listing  A.l:  SourceCode/BATCAALFE.m 

1  function  St _dot =BATCAM_FE ( St ) 

'/.this  calculates  the  first  order  derivatives  for  state  pr  opogat  i  on .  .  . 
given 

7, the  state  and  inputs  using  6  DOF  flat  earth  equations  (Lewis& 

7.  Stevens  p46)  The  state  (St)  is  a  1x17  state: 

7.  [u  v  w]  velocities  in  the  body  frame  (1st  three  elements) 

6  7.  [p  q  r]  angular  velocities  in  the  body  frame  (2nd  three  ... 

elements ) 

7.  [phi  th  psi]  attitude  (roll  pitch  yaw)  (3rd  three  elements) 

7.  [pe  pn  pd]  position  in  the  NED  frame  (4th  three  elements) 

7. 

7.  [Th  d_0  d_T  We  Wn] 

11  7.  Th  =  throttle  in  percent  e.g.(.5  =  507.) 

7.  d_0  opposing  elevon  cmd  in  degrees 

7.  d_T  tandem  elevon  cmd  in  degrees 
7.  We  is  east  component  of  wind 
7.  Wn  is  north  component  of  wind 

16 

7«separate  out  State  and  Inputs 
Inp  =  St (13: 17)  ; 
vb  =  St  (1:3)  ; 
wb=St (4:6) ; 

21  att  =  St  (7:9); 
pNED  =  St  (10  :  12)  ; 

Th=Inp (1)  ; 
d_0  =  Inp (2)  ; 
d_T  =  Inp (3)  ; 

26  We  =  Inp  (4)  ; 

Wn= Inp (5)  ; 

7«mass  properties 
m= . 425 ; 
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31  g=9 . 81 ; 

lxx=0 . 0028 ; 

Iyy =0 . 0088 ; 

Izz  =0 . 009  1  ; 

J=[Ixx  0  0 ; 0  Iyy  0;0  0  Izz]; 

36 

7,  calculate  necessary  Matrices 

0m_b=[0  -wb(3)  wb(2);wb(3)  0  -wb  ( 1 )  ; -wb  (2)  wb(l)  0]; 
Bb=[cos(att(2))*cos(att(3))  cos(att(2))*sin(att(3))  -sin(att(2)) ; 

-cos (att (1) ) *  sin ( att (3) )+sin(att (1) ) *  s in ( att (2) ) *cos (att (3) )  . 

cos (att (2) )*cos(att (3) )+sin(att (1) )*sin(att (2) )*sin(att (3) ) 
sin  (  att  (  1 )  )  *  cos  (  att  (2)  )  ;  ... 

41  sin(att(l))*sin(att(3))+cos(att(l))*sin(att(2))*cos(att(3)) 

sin ( att (1) ) *cos (att (3) )+cos (att (1) ) *sin (att (2) ) *  s  in (att (3) ) 
cos(att(l))*cos(att(2))] ; 

EP=[1  tan ( att (2) ) *  sin ( att ( 1 ) )  tan ( att (2) ) *  cos ( att ( 1) )  ;  ... 

0  cos (att  (1)  )  -  sin ( att  (  1 )  )  ;  ... 

0  sin  (  att  (  1 )  )  /  cos  (  att  (2)  )  cos  (  att  (  1 )  )  /  cos  (  att  ( 2)  )  ]  ; 

46  '/.calculate  airspeed  (assumes  0  down  component  of  wind) 

”/0vr  =  velocity  relative  to  wind 
vr=vb-Bb* [Wn ;  We;  0]; 

Va=norm ( vr ) ; 

'/.calculate  angle  of  attack,  side  slip  angle 
51  alph=atan (vr (3) / vr (1) ) ; 
beta=asin ( vr (2) /Va) ; 

7,  call  function  that  calculates  aero  forces 
if  Th  <  0 
56  Th=0 ; 

end 

if  Th  >  1 
Th  =  1  ; 

end 
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Fa=  AeroForces  (Va ,  alph  ,  beta,  Th  ,  d_0  ,  d_T); 

Fxyz  =  Fa  (1:3)  ’  ; 

LMN  =  Fa  (4:6)  1  ; 

66 

7,  Calculate  time  derivatives 
d_vb  =  -Om_b*vb  +  Bb*  [0  0  g]  ,+Fxyz/m; 
d_wb=-J~-l*Om_b*J*wb+J~-l*LMN  ; 
d_att=EP*wb ; 

71  d_pNED  =Bb 1 *vb ; 


St_dot= [d_vb ;  d_wb ;  d_att ;  d_pNED]; 

Listing  A. 2:  SourceCode/AeroForces.m 

function  Fa= AeroForces (Va , alpha ,  beta,  Th ,  d_0 ,  d_T) 

‘/.This  function  calculates  the  aerodynamic  forces  (Fxyz)  and  moments... 
( LMN  on 

‘/.the  BATCAM  airframe.  All  units  are  in  metric 
7,  Input  s 

5  7.  Va=airspeed 

7.  alpha=  angle  of  attack 

7.  beta=  side  slip  angle 

7.  Th  =  throttle  percentage 

7.  d_0  =  opposing  elevon  deflection  command  (in  degrees) 

10  7.  d_T  =  tandem  elevon  deflection  command  (in  degrees) 

7«  Out  puts 

7.  Fa  =  [Fx  Fy  Fz  L  M  N] 

7.  [Fx  Fy  Fz]  are  the  forces  in  the  xyz  direction 
7.  [L  M  N]  are  moments  about  the  xyz  axes 

15 

7o  BATCAM  properties 

rho  =  1.204;  7«air  density  in  metric 
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S  =  8 . 73*5 . 94*2*2 . 54~2/100~2  ;  ‘/.wing  area  [m~2] 
b  =  2 1  *  2 . 54/100 ; 

20  c_bar  =5 . 94*2 . 54/ 100  ;  '/.mean  geometric  chord  [m] 
q_bar=rho/2*Va~2 ;  '/.free  stream  dynamic  pressure 
d  =  8. 5*2. 54/100;  "/prop  diameter 


25  7, Calculate  coefficients 
CL=5 . 73* alpha  +  0.75; 

CD=1 . 203*alpha  +  0.075; 

CY  =  -0 . 573*  beta ; 

Clmn_T=  [0 . 000  .0167  -.00]; 

30  Clmn_0=  [  . 0021  .00  -.0072]; 

7,  Calculate  forces/moments 
D  =  q_bar  *S*CD  ;  /.Drag 
L  =  q_bar  *S*CL  ;  /.Lift 
35  Y  =  q_bar  *S*CY ;  7»Side  force 

LMN_T  =  q_bar  *S*  [b  c_bar  b]  .  *  Clmn_T  *  d_T  ;  7«Tandem  elevon  moments 
LMN_0  =  q_bar  *S*  [b  c_bar  b]  .  *Clmn_0*d_0  ;  7«0pposing  elevon  moments 

7«  Calculate  Thrust 
40  nu  =  117*Th;  7«revolutions/s 
J  =  Va/(nu*d);  7«advance  ratio 
Ct  =  -.1*J+.135;  7»  coefficient  of  thrust 
T  =  Ct*rho*nu~2*d~4;  "/.Thrust 

45  7. output 

Fa= [T-D  Y  -L  LMN_T  +  LMN_Q ]  ; 
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Appendix  B.  CC  Algorithm  Source  Code 

The  C++  source  code  for  the  dialog  box  that  implements  the  CC  algorithm  that 
interfaces  with  VC  is  on  the  accompanying  CD.  Only  the  Matlab®  files  that  imple¬ 
ment  the  algorithm  are  included.  The  C++  code  utilizes  the  concepts  in  Procerus’s 
DevDemo  sample  application  which  illustrates  how  to  construct  and  transmit  VC 
packets,  and  also  the  Matlab®  capability  that  interfaces  with  C++  using  the  “mcc 
-B  csharedlib:***  ***.m  ...”  command. 

B.l  Matlab  Files  for  CC  Algorithm 

Here  are  the  associated  Matlab®  routines  and  dependencies  implemented  in  the 
C++  source  code.  Refer  to  Section  3.6.4  for  descriptions  of  the  routines. 

Listing  B.l:  SourceCode/EnterOrbit.m 

function  [LatLonHdg  r] =EnterOrbit (  LatLon_t  ,  alt,  th_d  ,  Vg  ,  numUAV ,  offst, 
GPS_ht ) 

'/.This  Function  provides  a  matrix  (nx3)  of  the  [lat  Ion  hdg]  for  n  UAV  ’  s  to 
'/.enter  an  orbit  around  a  target  positioned  at  LatLon_t  =  [lat  Ion]  .  The 
4  '/.orbit  is  a  circle  divided  into  numUAV  equal  parts  at  altitude  (alt)  and 
"/.with  sensor  depression  angle  th_d  .  All  vehicles  enter  orbit 
7,  counter  clockwise  . 

7,  Input  s  : 

1  LatLon_t  =  [lat  Ion]  in  degrees  for  the  target  position 
9  /  alt  =  altitude  of  orbit  (in  meters  above  ground  level) 

7.  th_d  =  sensor  depression  angle  (radians).  Algorithm  assumes  the  sensor 
7.  is  side  facing  to  view  the  target  in  the  middle  of  the  orbit 

7o  Vg  =  the  ground  speed  of  the  vehicle  (user  must  find  ground  speed 

7.  given  airspeed  and  wind  conditions) 

14  7.  numUAV  =  the  number  of  vehicles  to  orbit  the  target 

7«  offst  =  any  angular  offset  (in  radians).  At  offst=0,  the  first  point 
7.  is  at  0  rad  (at  point  directly  east  of  target  location) 

1  GPS_ht  =  height  above  mean  sea  level  (in  meters) 

7o  Out  puts  : 

19  7.  LatLonHdg  =  nx3  matrix  where  each  row  denotes  the  [lat  Ion  hdg]  to 
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7,  enter  the  orbit.  n  equals  numUAVs.  lat/lon  are  in  decimal  degrees, 

7.  heading  is  in  radians  . 

7,  r  =  orbit  radius  (in  meters) 

7. 

24  7.  This  function  needs  the  Orbit_r  function  ,  Cart2GPS  function,  and 
7.  Zeroto2pi  function. 

[r  phi]=  Orbit_r(alt,  th_d  ,  Vg); 
ii=l : numUAV ; 

29  theta  =  2*pi/numUAV*  ( ii -1)  ;  7«create  angle  increments 

xy  =  r*  [cos  (theta  +  of  f  st  )  ;  sin  (  thet  a  + of  f  st  )  ]  ;  7«create  xy  coordinates 
[lat  Ion] =Cart2GPS (LatLon_t  ,  xy(l,  :)  ’,  xy(2,:)’,  GPS_ht); 
hdg  =  Zeroto2pi (2*pi-(theta+offst) )  1  ; 

LatLonHdg = [ lat  Ion  hdg] ; 

Listing  B.2:  SourceCode/CreateOrbit.m 

function  [Wp  d_mtx  Hdg_mtx  rf ] =CreateOrbit (  GPS,  Hdg,  r_min  ,  GPS_ht) 

7»This  function  creates  a  series  of  waypoints  based  on  initial  and  final  GPS 
70points  and  headings.  The  lengths  of  all  paths  are  the  same  within  5  meters. 
5  70The  waypoints  are  along  a  Dubins  path.  A  waypoint 

7ois  provided  at  least  every  90  degrees  along  curved  portions  of  the  path 

7. 

7,  Input  s  ar  e  : 

7o  GPS  =  nx4  matrix  denotes  the  initial  and  final  values  of  GPS 
10  7.  coordinates  (lat,  Ion)  for  each  vehicle  (in  degrees).  Each  row  is  the 
7.  set  for  each  vehicle  [Lat_initial  Lon_initial  Lat_final  Lon_final] 

7.  Hdg  =  nx2  matrix  denoting  the  initial  and  final  headings  for  each 
/  vehicle  (in  radians).  Each  row  is  for  each  vehicle  [Hdg_i  Hdg_f] 

7.  r  =  the  minimum  turn  radius  (scalar--in  meters) 

15  7«0utputs  are: 

7»  Wp  =  2nxm  matrix  providing  the  waypoints  (in  degrees).  Rows  1  and  2  ... 
denote  the  set 

7o  of  points  for  vehicle  l[Lat;  Lon],  Rows  3  and  4  denote  the  set  for  ... 
vehicle  2 ,  and 
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25 

30 

35 
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45 

50 


7.  so  on.  Since  the  number  of  columns  for  each  vehicle  can  vary,  the 
”/o  dimension  m  denotes  the  length  of  the  longest  set.  The  end  of  shorter  . 
sets  are 

7.  padded  with  the  number  -100  since  latitude  cannot  exceed  +/-90.  The 
7.  user  must  account  for  the  fact  that  if  he /  she  sees  a  -100  in  a  latitude 
7o  measurement  ,  the  end  of  the  list  has  been  reached. 

7.  d  =  nxm  matrix  denoting  the  distance  at  each  waypoint  of  the  Dubins 

l  path(in  meters)  of  each  vehicle -total  dist  is  at  last  point.  this  is 

7»  also  padded  with  -100  for  shorter  vectors 

i  Hdg_Mtx  nxm  matrix  of  heading  at  each  waypoint 

7.  r  =  lxn  vector  denoting  the  radius  of  the  Dubins  path.  Realize  this  . 
algorithm  assumes 

i  that  the  radii  for  each  end  of  the  path  are  the  same. 


ro=r_min ; 

numUAV=size (GPS) ;  numUAV=numUAV (1) ; 

7.  i=3020  :  40  :  3160  ; 

7.  GPSi  =  r _GPS ( i  ,  :  )  ; 

i  [lat_f  lon_f]  =  Cart2GPS(  [39.3435  -86.0345],[r_orbit  0  -r_orbit  0]  ’  ,  .  .  . 
7.  [0  r_orbit  0  -r_orbit]  ’  ,216)  ; 

/  GPSf  =  [ lat_f  lon_f ]  ; 

7.  Hi  =  att  (i  ,  3)  ; 

7.  Hf  =  [0  270  180  90]  ’*pi/180; 

/  GPS_ht=Alt_GPS (1) ; 

7.  ro  =20  ; 

rf=zeros (1 , numUAV ) +r _min ; 


7«  generate  initial  distances 
for  ii=l : numUAV 

[W  d_temp  Htmp]=Dubins([-l  -1 ; 0  0],[GPS(ii,l:2);GPS(ii,3:4)],... 

[Hdg ( ii  ,  1 )  Hdg ( ii  , 2) ]  ,  ro  ,  GPS_ht); 
d ( i i ) =d_t emp ( 1 ) ; 
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if  ii  >  1  &&  d ( i i ) >dmax ( 1 ) 


dmax=d_temp ; 
Wmax=W ; 
Hmax=Htmp ; 
elseif  ii  ==  1 

dmax=d_temp ; 
Wmax=W ; 
Hmax=Htmp ; 

end 


end 

'/.sort  from  largest  to  smallest 

[ds  iis]=sort(d, ’descend') ; 

lenWmax  =  size (Wmax)  ;  lenWmax  =  lenWmax  (2)  ; 

Wp=zeros (2*numUAV , lenWmax ) ;  %  create  initial  Wp  matrix  with  zeros 

Wp  (  (2* i is ( 1 ) - 1 )  :(2*iis(l))  ,:)=Wmax; 

d_mtx  =  zeros  (numUAV  ,  lenWmax);  '/.create  initial  d_mtx  w/zeros 

d_mtx(iis(l)  ,:)=dmax;  ’/.place  max  dist  row  in  right  place  according  to  sort 

Hdg_mtx  =  zeros  (numUAV  ,  lenWmax);  "/.create  initial  Hdg_mtx  w/zeros 

Hdg_mt  x  ( i  is  ( 1 )  ,:)=Hmax;  '/.place  max  dist  row  in  right  place  according  to  sort 


Kp= . 05 ; 


"/.adjust  radii  to  get  distances  all  the  same 
for  ii=2: numUAV 
ctr  =  1 ; 

dd  =  100;  "/.get  things  started 
dr  =2  ; 


d_dr  =  2 ; 
extra_turn=0 ; 
extra_turn_flg=l ; 
dd_sign=l ; 
ctr_f lg=0 ; 
while  abs(dd)  >  5 

rf (iis (ii) ) =rf (iis (1) ) +dr  ; 

7. recalculate  with  new  radius 
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[W  d_t emp  Hdg_tmp]=Dubins([-l  -1 ; extra.turn  0],[GPS(iis(ii),l:2);GPS(... 
iis(ii)  ,3:4)] 

[Hdg (iis(ii)  ,1)  Hdg ( iis ( ii  )  ,  2)  ]  ,  rf ( iis ( ii )  )  ,  GPS_ht); 
d  (  iis  (  ii  )  )  =d_temp  (  1 )  ;  7, grab  path  length 

if  extra.turn  &&  extr a_turn_f lg 
extra_turn_flg=0; 
if  d(iis(ii))  >  dmax 
d_dr=-d_dr ; 

end 

end 

dd  =  d  (  iis  ( 1 )  ) -d  (  iis  (  ii  )  )  ;  70find  difference  btwn  longest  and  current 

ctr=ctr +1 ; 
dr=dr +d_dr ; 

if  (dd*dd_sign  <  0)  7.  detects  an  overshoot 

dd_sign=-dd_sign ; 

d_dr  =  -d_dr /2 ;  7ohalve  increment  of  change 

end 

if  (ctr  >  100)  7.  if  doesn’t  converge  add  extra  turn  and  reset  ... 

everything 

if  ctr_flg  7«if  it  doesn’t  converge  after  200  just  return 
return ; 

end 

extra_turn=l ; 
dr  =0  ; 
d_dr =1 ; 
ctr  =0 ; 
ctr_f lg=l ; 
dd_sign=l ; 

end 


end 

lenW=size (W) ;  lenW=lenW (2) ; 
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if  lenW  >  lenWmax 

120  Wp=[Wp  (  zeros  ( 2*  numUAV  ,(  lenW  -  lenWmax  ))-  100 )]  ;  "/.add  length,  to 

d_mtx=[d_mtx  zeros (numUAV ,  lenW - lenWmax )- 100] ;  "/.add  length  to 
Hdg_mtx=  [Hdg_mtx  zeros  (numUAV  ,  lenW  -  lenWmax  )- 100]  ;  */,add  length 
Hdg_mtx 
lenWmax=lenW ; 
elseif  lenW  <  lenWmax 

125  W=[W  (zeros  (2  ,  (  lenWmax  -  lenW ))  — 100)]  ;  "/.extend  W 

d_temp=  [d_temp  zeros  (  1  ,  lenWmax  -  lenW)  -  100]  ;  '/.extend  d_temp 
Hdg_tmp=  [Hdg_tmp  zer  os  (  1  ,  lenWmax  -  lenW  )- 100]  ;  "/.extend  Hdg_tmp 

end 

Wp((2*iis(ii)-l):(2*iis(ii)),:)=W;  7.  insert  W  into  Wp 

130  d_mtx  (  iis  (  ii  )  ,  :  )  =d_temp  ;  "/.insert  d_temp  into  d_mtx 

Hdg_mtx(iis(ii)  ,  :) =Hdg_tmp ; 


Wp 

d_mtx 
to  . 


end 


Listing  B.3:  SourceCode/UpdateGrndSpd.m 

function  Vd=Updat eGrndSpd ( Cur_GPS ,  Cur_Wp ,  Wp_mtx ,  d_mtx ,  V_nom ,  GPS_ht) 

7.  This  function  calculates  the  necessary  ground  speed  to  complete  the 
7.  remaining  waypoints  in  the  flightplan  so  that  all  vehicles  arrive  at  the 
4  7,  last  waypoint  at  the  same  time. 

7.  Input  s  : 

7.  Cur_GPS  =  nx2  array  of  [lat  Ion]  coordinates  of  each  Vehicle.  n 
7,  denotes  the  number  of  vehicles. 

7.  Cur_Wp  =  lxn  vector  of  numbers.  Each  number  denotes  the  next  waypoint 
9  7,  in  the  flight  plan. 

7.  Wp_mtx  =  2nxm  matrix  where  every  row  pair  are  the  waypoints  for  each 

7.  vehicle.  See  the  CreateOrbit  command  for  more  details. 

7,  d_mtx  =  nxm  matrix  of  distances  remaining  in  the  flightplan.  See  the 
7.  CreateOrbit  command  for  more  details 
14  7.  SpdRnd  =  1x2  vector  that  contains  the  min  and  max  velocity  [Vmin  Vmax] 

7,  GPS_ht  =  the  reference  GPS  height  ( height  above  MSL) 

7. Output s  : 
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7,  Vd  =  lxn  vector  of  updated  velocities  [m/s]  required  for  each  vehicle  to 
70  arrive  all  at  the  same  time.  The  routine  uses  the  Vmax  value  for  the  ... 
flight 

19  7,  plan  with  the  greatest  distance.  If  no  solution  exists  within  the  given 
1  velocity  range,  the  function  returns  a  Vd  vector  of  zeros. 

7,  Calculate  distances  from  current  position  to  next  waypoint 
numUAV= length ( Cur _Wp ) ; 

24  numWp  =  s ize ( Wp_mtx )  ;  numWp  =  numWp  (2)  ; 
for  ii=l:numUAV 

[d_tmp  H]=DistFromGPS([Cur_GPS(ii,:);  Wp_mtx((2*ii-1)  : (2* ii )  , Cur_Wp ( ii ) ) . .  . 
’] , GPS_ht ) ; 

d_tot ( ii) =d_tmp (2) +d_mtx (ii , Cur_Wp (ii) ) ; 

end 

29  d_mean=mean ( d_tot ) ; 

t  =  d_mean  /  V_nom ;  7«find  the  time  to  complete  using  nominal  V  and  mean  dist. 

Vd=  d_t ot / 1  ; 


Listing  B.4:  SourceCode/Airspeed.m 

function  Va= Air speed (Vg ,  Vw ,  V_limits) 

1  This  function  returns  the  necessary  airspeed,  necessary  heading,  crab 
3  7.  angle,  and  a  flag  stating  whether  the  needed  airspeed  exceeds  limits. 

7.  Input  s  : 

1  Vg  =  1x2  vector  that  contains  desired  [Ground_Speed  Heading (rad) ] 

7«  Vw  =  1x2  vector  that  contains  [Wind_Speed  Wind_Heading  (rad)  ] 

7.  V_limits  =  1x2  vector  with  max/min  values  [V_max  V_min] 

8  7« Outputs 

7o  Va  =  1x4  vector  that  contains  [Airspeed,  Heading  (rad)  ,  crab.angle  (rad)  , 
°l  flag] 

7.  All  speeds  must  be  in  the  same  units,  all  angles  must  be  in  radians. 

7«  Positive  crab  angles  denote  CCW  angle  from  ground  heading  to  airspeed 
13  /  heading . 

th_g=pi/2-Vg (2) ; 
th_w=pi/2-Vw (2) ; 

Vg_vec=Vg(l)*[cos(th_g)  sin(th_g)] ; 
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Vw_vec=Vw(l)*[cos(th_w)  sin(th_w)] ; 

18  Va_vec =Vg_vec -Vw_vec ; 

Va(l)=norm(Va_vec)  ;  “/.airspeed 

Va(2)=Zeroto2pi(pi/2-atan2(Va_vec(2) ,  Va_vec(l))) ; 
Va(3)=atan2(Va_vec (2) ,  Va_vec (1) ) -th_g ; 
if  (Va(l)  >  V_limits(l))  I  I  (Va(l)  <  V_limits(2)) 
23  Va(4)=l;  “/.calculated  airspeed  exceeds  limits 

else 

Va(4)=0;  “/.Calculated  airspeed  is  good 

end 


Listing  B.5:  SourceCode/MaintainOrbit.m 

1  function  Vd=MaintainOrbit (GPS ,  GPS_tgt , R_Vo ,  Kp_Slop ,  GPS_ht) 

7, This  function  is  used  to  maintain  a  CCW  orbit  around  a  specified  target 
“/.location.  The  output  is  a  set  of  ground  speeds  that  will  adjust  the 
“/.current  positions  of  the  vehicles  to  maintain  an  equally  spaced  orbit. 

7,  Input  s  : 

6  7.  GPS  =  nx2  array  of  current  lat  itude  /  longitude  positions  of  each 

7,  vehicle.  Each  row  is  [lat  Ion]  and  n  is  the  number  of  vehicles.  All 
7,  lat  itude  /  longitude  values  are  in  decimal  degrees. 

7.  GPS_tgt  =  1x2  vector  of  the  target  location,  [lat  Ion] 

7.  R_Vo  =  1x2  vector  of  [Orbit_Radius  Nominal_Speed]  Radius  is  in  meters, 

11  7,  Speed  is  in  m/s  . 

7.  Kp_Slop  =  1x2  vector  of  [Proportional.Gain  Angular_Slop]  .  Kp  is  a 

7,  mult  iplier  on  the  angular  error  ( e  =  Ang  -  Ang_des)  Slop  is  the  angular 

7,  window  that  is  considered  in  the  correct  position,  (e.g.  0.1745  rad  (10 

7.  deg)  is  +/-  5  degrees  is  considered  ’good')  This  is  in  radians. 

16  7. Outputs  : 

7.  Vd  =  lxn  vector  of  desired  velocities  (in  m/s) 

7,  Calc  distances  from  center 

numUAV  =  size (GPS)  ;  numUAV  =  numUAV  (1)  ; 

21  d  =  DistFromGPS  (  [GPS_tgt  ;  GPS],  GPS_ht);  “/.Calculate  distances 


d=d (2 : numUAV ) ; 

7,  Calc  angles  from  center 
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[pn  pe ] =  GPS2Cart ( [GPS_tgt  (1)  ; GPS ( :  , 1) ]  1  ,  [GPS.tgt  (2)  ; GPS  (  :  ,2) ]  ’  , GPS_ht)  ; 


26  UAV_ang  =  atan2  (pn  ,  pe  )  ;  “/.calc  angular  position  for  each.  UAV 
UAV_ang=UAV_ang (2 : (numUAV+1) ) ; 

dUAV_ang  =  UA¥_ang -UAV_ang  ( 1 )  ;  “/.adjust  WRT  first  UAV  os  UAV1  is  at  0  rad 
dUAV_ang=Zeroto2pi (dUAV_ang) ; 

31  [dUAV_ang_s  iis  ]  =  sort  ( dUAV_ang  ,  ’ascend’);  7.  sort  angle  differences 
ii=l : numUAV ; 

Ang_Des  =  (  ii  -  1 )  *2* pi /numUAV  ;  “/.calculate  desired  angular  position 
Ang_Err = Ang_Des - dUAV_ang ( i is ) ; 

Vd=zeros (1 , numUAV )+R_Vo(2) ; 

36  for  i  i  =  2 :  numUAV  "/.loop  thru  sorted  list  and 
if  abs ( Ang_Err ( iis  (  ii  )  )  )  <  Kp_Slop (2) 

Vd(iis(ii))=R_Vo(2) ; 

else 

Vd(iis (ii) ) =R_Vo (2) +Kp_Slop (1) *Ang_Err (iis (ii) ) ; 

41  end 

end 


Listing  B.6:  SourceCode/Dubins.m 

7.7.  Dubins  path  calculator 

2  function  [Wp  d  Hdg_W] =Dubins (turn ,  GPS,  Hdg ,  r,  GPS_ht) 

7.  Dubins  path  calculator 

7.  This  function  calculates  a  set  of  GPS  waypoints  and  total  distance  for 
7.  a  turn  -  straight  -  turn  path(Dubins  path)  given  an  initial  position  and 
7.  heading  . 

7  7.  Inputs  : 

7.  turn  is  a  2x2  vector  denoting  direction  of  turns.  The  first  row 

7.  designates  direction  of  turns,  first  number  of  first  row  denotes 

7.  initial  turn,  second  denotes  final  turn. 

7.  For  the  first  row 

12  7.  [1  1]  =  CCW  CCW  turn  scheme  ,  [-1  -1]  =  CW  CW  turn  scheme  ,  no 

7.  extra  turns 

7.  [-1  1]  =  CW  CCW  turns,  [1  -1]  =  CCW,  CW  turns 
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°i  The  second  row  denotes  whether  to  add  an  extra  turn  to  either  end 

1  [0  0]  denotes  no  extra  turns 

17  °/o  [1  0]  denotes  extra  turn  on  initial  circle 

7,  [0  1]  denotes  extra  turn  on  final  circle 

70  total  turn  mtx  for  CCW  CCW/extra  no  extra  is  [1  1 ;  1  0] 

%  GPS  is  a  2x2  matrix  of  GPS  coordinates  of  initial  and  final 
°i  coordinates  . 

22  1  [Lat_initial  Lon_initial] 

7,  [Lat_final  Lon_final  ] 

7o  Hdg  is  a  1x2  vector  with  initial  and  final  headings  in  radians 
7.  [Heading_  initial  Heading_final] 

7.  r  is  the  radius  of  the  turns  .  r  and  d  are  in  meters 
27  7.  GPS_ht  is  the  reference  height 

7. 

°l  Output s  : 

7,  Wp  is  a  2xn  set  of  GPS  waypoints.  One  waypoint  is  provided  for  every 
7.  90  deg  of  circle  sweeped  ,  one  at  the  final  circle  entry,  one  at  the 

32  7.  final  destination,  and  one  past  the  destination  in  the  direction  of 

/  the  final  heading????? 

7.  d  lxn  vector  that  denotes  the  distance  remaining  at  each  waypoint. 

7,  d(l)  is  the  total  distance,  d(n)  is  zero. 

7,  A  is  the  total  radians  swept  in  the  path. 

37 

7,  convert  GPS  to  meters 

[y  x] =GPS2Cart (GPS ( : , 1) ,  GPS(:,2),  GPS_ht); 

7,  convert  Heading  to  NED  coordinates  (  +  x  =  East,  +y  =  North)  where  zero 
42  7.  angle  is  along  x  axis 
NEDHdg=pi/2-Hdg ; 

7.  find  turn  circle  centers 
7o  c  i  r  c  1  e  1 

47  x_cl=x(l)+r*cos(NEDHdg(l)+turn(l  ,  l)*pi/2)  ; 

y_cl=y(l)+r*sin(NEDHdg(l)+turn(l  ,  l)*pi/2)  ; 

7o  c  i  r  c  1  e  2 
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x_c2=x(2)+r*cos(NEDHdg(2)+turn(l  ,2)*pi/2)  ; 
y_c2=y(2)+r*sin(NEDHdg(2)+turn(l  ,2)*pi/2)  ; 

52 

7.  find  points  of  tangency  between  two  circles 

if  turn  ( 1  ,  1 )  *  turn  ( 1  , 2)  ==  1  “/.for  initial  and  final  circles  in  the  same  dir 

if  turn(l,l)  ==  1  70for  ccw  on  circle  1  and  circle  2 
Ai_cl=Zeroto2pi (NEDHdg (1) +3*pi/2) ; 

57  Af _c 1 = Zer ot o2pi ( at an2 (y_c2 -y_c 1 ,  x_c2 -x_c 1 ) +3* pi /2 ) ; 

Ai_c2=Af _cl ; 

Af_c2=Zeroto2pi (NEDHdg (2) +3*pi/2) ; 

Al=Zeroto2pi (Af_cl -Ai_cl ) +turn (2 , 1) *2*pi ; 

A2=Zeroto2pi (Af_c2-Ai_c2)+turn(2 ,2) *2*pi ; 

62  else  70  cw  for  circle  1  and  cw  for  circle  2 

Ai_cl=Zeroto2pi (NEDHdg (l)+pi/2) ; 

Af_cl=Zeroto2pi (at an 2 (y_c2-y_cl ,  x_c2-x_cl)+pi/2) ; 

Ai_c2=Af _cl ; 

Af_c2=Zeroto2pi (NEDHdg (2) +pi/2) ; 

67  Al=-m2piToZero(Af_cl-Ai_cl)+turn(2,l)*2*pi; 

A2=-m2piToZero(Af_c2-Ai_c2)+turn(2,2)*2*pi; 

end 

d_str =norm ( [ x_c2 - x_c 1  y_c2-y_cl]) ; 
else  7«for  circles  in  opposite  direction 
72  7«check  if  points  are  too  close  for  this 

if  norm ( [x_c2 -x_cl  y_c2-y_cl])  <  2*r 

sprintf ( 1 \nPoints  are  too  close  for  cw/ccw  turn  —  abort .  ; 

return 

end 

77  if  turn(l,l)  ==  1  70  for  ccw  for  circle  1,  cw  for  circle  2 

7oCompute  angle  btwn  circle  centers  and  tangent  line 
th=atan(r/norm( [x_c2-x_cl  y_c2-y_cl] ) ) ; 

Ai_cl=Zeroto2pi (NEDHdg (1) +3*pi/2) ; 

Af_cl=Zeroto2pi(atan2(y_c2-y_cl ,  x_c2-x_cl)+3*pi/2+th) ; 

82  Ai_c2=Zeroto2pi ( atan2 (y_c2 -y_cl ,  x_c2 -x_c 1 ) +pi /2+ th ) ; 

Af_c2  =  Zeroto2pi (NEDHdg (2) +pi/2)  ; 

Al=Zeroto2pi (Af_cl -Ai_cl ) +turn (2 , 1) *2*pi ; 
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A2=-m2piToZero (Af_c2-Ai_c2)+turn(2 ,2) *2*  pi ; 
else  '/.find  angles  for  cw  circle  1  and  ccw  circle  2 
87  th=atan (r/norm ( [x_c2 -x_cl  y_c2 -y_c 1 ] ) ) ; 

Ai_cl=Zeroto2pi (NEDHdg(l)+pi/2) ; 

Af_cl=Zeroto2pi (atan2 (y_c2-y_cl ,  x_c2-x_cl)+pi/2-th) ; 
Ai_c2  =  Zeroto2pi(atan2(y_c2-y_cl  ,  x_c2-x_cl)+3*pi/2-th) 
Af_c2=Zeroto2pi (NEDHdg (2) +3*pi/2) ; 

92  Al=-m2piToZero(Af_cl-Ai_cl)+turn(2,l)*2*pi; 

A2=Zeroto2pi (Af_c2-Ai_c2)+turn(2 ,2) *2*pi ; 

end 

d_str=norm([(x_c2+r*cos(Ai_c2))-(x_cl+r*cos(Af_cl))  .  .  . 
(y_c2+r*sin(Ai_c2) ) -(y_cl+r*sin(Af_cl) )] ) ; 

97 

end 

7,  Calculate  A 
A=A1 +A2 ; 

102 

7,  Calculate  waypoint  s 

7,  calculate  points  on  circle  1 
Aldiv  =  ceil ( A 1 / (pi/2) )  ; 
ii = 1 : (Aldiv+1) ; 

107  x_wp  =  x_c 1 +r *  cos ( Ai_c 1 + A1 *  turn ( 1 , 1 ) * ( ii - 1 ) / A ldi v )  ; 

y_wp=y_cl+r*sin(Ai_cl+Al*turn (1 , 1) *  (  ii -1) /Aldiv) ; 

d=r*Al*(ii-l)/Aldiv; 

if  turn (1,1) ==1 

Hdg_W=Zeroto2pi (2*pi-(Ai_cl+Al*(ii-l) /Aldiv) ) ; 

112  else 

Hdg_W=Zeroto2pi (pi-(Ai_cl -Al*(ii -1) / Aldiv) ) ; 

end 

/•calculate  points  for  circle  2  (including  endpoint 
117  dii=length (x_wp) ; 

A2div  =  ceil (A2/ (pi/2) )  ; 
ii = 1 : (A2div+1) ; 
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x_wp=[x_wp  x_c2+r*cos(Ai_c2+A2*turn(l,2)*(ii-l)/A2div)] ; 
y_wp=[y_wp  y_c2+r*sin(Ai_c2+A2*turn(l,2)*(ii-l)/A2div)] ; 

122  d(ii+dii)=r*Al+d_str+r*A2*(ii -1) / A2div ; 

if  turn ( 1 , 2) ==1 

Hdg_W(ii+dii)=Zeroto2pi (2*pi-(Ai_c2+A2*(ii-l)/A2div)) ; 

else 

Hdg_W (ii+dii)=Zeroto2pi (pi-(Ai_c2-A2*(ii-l) / A2div) ) ; 

127  end 

"/convert  to  GPS  coordinates 
d=d(length(d)  )  -  d; 

132  [Wp ( 1  ,  : )  Wp (2 , : ) ] =Cart2GPS ( [GPS (1  ,  1)  GPS(1,2)],  x_wp ,  y_wp ,  GPS_ht); 

Listing  B.7:  SourceCode/Orbit _r.m 

function  [r  phi]  =  Drbit_r (h ,  Theta_s,  Vg) 

2  '/.This  function  solves  for  the  orbit  radius  and  bank  angle  of  an  orbiting 
70airplane/UAV  given  an  altitude  (h),  a  side  sensor  depression  angle 
°/(theta_s),  and  a  ground  speed(Vg).  This  function  assumes  that  the  bank 
"/angle  is  small  and  only  returns  the  larger  radius  (the  solution  actually 
'/.has  two  solutions)  All  answers  are  given  in  meters  and  radians.  Solution 
7  "/also  assumes  that  altitude  is  constant. 

"/  Input  s  : 

"Z  h  =  altitude  [meters] 

"/  Theta_s  =  depression  angle  of  side  sensor 
"/  Vg  =  ground  speed 
12  %  Output  s  : 

"/  r  =  orbit  radius 
"Z  phi  =  bank  angle 

"/solve  the  quadratic  equat  ion --give  only  larger  solution 
17  g  =  9.81;  "/  acceleration  of  gravity 

a=  1 ; 

b  =  -(h*tan ( Thet  a_  s )+Vg~2/g*tan(Theta_s))  ; 
c  =  -h*Vg  ~2/g ; 
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r=(-b+sqrt (b~2-4*a*c))/(2*a) ; 

22  phi=atan(Vg~2/ (g*r) ) ; 

Listing  B.8:  SourceCode/Cart2GPS.m 

7.  Cartesian  to  GPS  conversion 

function  [Lat  Lon] =Cart2GPS (GPS_o ,  x,  y,  GPS_ht) 

3  7,  This  function  generates  GPS  coordinates  [Lat  Lon]  from  a  given  GPS  origin 
7.  GPS_o  [Lat  Lon],  x  y  coordinates  (in  meters),  and  a  GPS  reference  height. 
7.  The  GPS  origin  must  correlate  with  the  xy  origin. 

7.  This  function  is  only  good  for  small  distances  (>1000km) 

8  7.set  the  origin 

Lat_o=GPS_o (l)*pi/180; 

Lon_o=GPS_o (2) *pi / 180 ; 

7.  Calculate  Rm  ,  Rp 

13  a  =  6378135;  7.equatorial  radius  in  meters 
e  =0 . 0818191908426 ; 

Rm=(a*(l-e~2) ) / (l-e~2*sin(Lat_o) ~  2)  .  '  1 . 5; 

Rp=a/ (l-e~2*sin(Lat_o) ~  2 ) ~  0 .5; 

18  Lat=Lat_o+y/(Rm+GPS_ht ) ; 

Lon=Lon_o+x/( (Rp+GPS_ht) *  cos (Lat_o) ) ; 

7o  convert  back  to  deg 
Lat=Lat*180/pi ; 

23  Lon  =  Lon *  180/pi ; 


Listing  B.9:  SourceCode/GPS2Cart.m 


7.7.  GPS2Cart 

2 

function  [pn  pe ] = GPS2Cart ( Lat ,  Lon,  h_o) 

7.  This  function  converts  a  set  of  latitude,  longitude,  height  values  to 
7.  cartesian  pn  ,  pe  values  in  meters.  The  first  point  is  set  as  the 
7.  origin.  Equations  come  from  slide  1-29  of  EENG  533  from  Spring  2008. 
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7  7,  this  assumes  Lat/Lon  are  in  degrees.  h_o  is  the  GPS  height  of  the  first 
°/„  point.  This  function  is  only  good  for  small  distances  (>1000km) 

"/.Convert  to  radians 
Lat=Lat*pi/180; 

12  Lon=Lon*pi / 180 ; 

7,  set  the  origin 
Lat_o=Lat (1) ; 

Lon_o  =  Lon  (  1 )  ; 


7, Calculate  Rm  ,  Rp 

a  =  6378135;  "/.equatorial  radius  in  meters 
e  =0 . 0818191908426 ; 

22  Rm=(a*(l-e~2))/(l-e~2*sin(Lat_o)~2)  .  ~  1 . 5  ; 
Rp=a/ (l-e~2*sin(Lat_o) ~  2) ~  0 .5; 
pe=(Rp+h_o) *cos (Lat_o) *(Lon-Lon_o) ; 
pn=(Rm+h_o)*(Lat-Lat_o) ; 
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Appendix  C.  Simulink  Model  Matlab®  Code 

This  appendix  contains  the  Matlab®  code  used  for  the  Simulink  model. 


C.l  BATCAM  Closed  Loop  Model 


Listing  C.l:  SourceCode/BAT CAM  S.rn 

function  dy=BATCAM_S (y ) 

”/o Calculates  the  derivatives  of  the  closed-loop  BATCAM  model  state 
%The  state : 

“/o  y(l)=  pN  North  position 

7«  y(2)=  pE  East  position 

7o  y(3)=  h  Altitude 

“/«  y  (4)  =  h  dot 

7«  y  ( 5 )  =  Va  Airspeed 

7o  y(6)=  R  Roll  angle 

“/«  y(7)=  P  Pitch  angle 

7«  y(8)=  Y  Yaw  Angle  (clock  angle) 

7.  Input  s 

“/«  y(9)=  Vac  Commanded  airspeed 

7«  y(10)=  Rc  Commanded  roll  angle 

7.  y  (  1 1 )  =  Pc  Commanded  Pitch  angle 
7.  y(12)=  wN  North  component  of  wind 

7«  y(13)=  wE  East  component  of  wind 

7,  Output  s 

7«  dy  (  1 )  =  pN  dot 

7.  dy  (2)  =  pE  dot 

7.  dy  (3)  =  h  dot 

7«  dy(4)=  h  dotdot 

7.  dy  (5)  =  Va  dot 

7.  dy  (6)  =  R  dot 

7«  dy  (7)  =  P  dot 

7.  dy(8)=  Y  dot 


7«First  Order  Constants 
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30  kV  =  1.3;  •/„  95"/0  of  commanded  in  2.3  seconds 

kR=2 . 3 ; 
kP= . 865 ; 

7.  roll  wind  disturbance 

35  kwR  =  l;  7.  1  =  1  degree  of  disturbance  for  every  m/s  of  wind 
wR  =  kwR*0.0174*norm([y(ll)  y  ( 12) ] )  ; 

7oBATCAM  properties 
m=0 . 425 ; 

40  g=9 . 81 ; 

rho  =  1.204;  7«air  density  in  metric 

S  =  8 . 73*  5 . 94*  2*2 . 54  ~  2/  1 00  “  2  ;  ‘/0wing  area  [m~2] 

7« b  =  2 1  * 2 . 54/100  ; 

70c_bar  =  5 . 94*2 . 54/100  ;  '/,mean  geometric  chord  [m] 

45  q_bar=rho/2*y  (5)  ~2 ;  7, free  stream  dynamic  pressure 
7od  =  8 . 5*2 . 54/100  ;  7«prop  diameter 

7.  Calculate  lift 
alpha=0 ; 

50  CL=5 . 73* alpha  +  0.7; 

L  =  q_bar *S*CL  ;  /.Lift 


dy(l)=y(5)*cos(y(8))*cos(y(7))+y(12); 

55  dy(2)=y(5)*sin(y(8))*cos(y(7))+y(13)  ; 
dy(3)=y(5)*sin(y(7)  )  ; 
dy  (4)  =L*cos  (y  (6)  )  *cos  (y  (7)  ) /m  -  g; 
dy(5)=kV*(y(9)-y(5)  )  ; 
dy (6)=kR*(y(10) -y (6) )+wR; 

60  dy (7)=kP*(y(ll) -y  (7) )  ; 
dy (8)=g/y(5) *tan(y(6) ) ; 

C.2  Mode  2  Controller 

Listing  C.2;  SourceCode/CC_Controllerl.m 
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function  Cmd=CC_Controller 1 (yin) 

'/.This  function  is  used  in  conj  unct  ion  with  Simul  ink  model  BATCAM_1  to  model 
3  '/.the  velocity  feedback  algorithm 
7,  Input  s  : 

1  yin(l:2)  =  wind  vector  [wNorth  wEast] 

%  yin(3:5)  =  position  of  UAV1  [pN  pE  Yaw] 

%  yin(6:8)  =  position  of  UAV2  [pN  pE  Yaw] 

8  l  yin(9:ll)  =  position  of  UAV3  [pN  pE  Yaw] 

1  yin(12:14)  =  position  of  UAV4  [pN  pE  Yaw] 

7o  Out  puts 

7.  [pNd  pEd]  =  desired  location  for  the  UAV 
7.  Chi_s  =  desired  path  heading 
13  7.  Vac  =  commanded  airspeed 
7.  he  =  commanded  altitude 

7.  Cmd  (1:5)  =  [pNd  pEd  Chi_s  Vac  he]  for  UAV1 
l  Cmd  (6:  10)  =  [pNd  pEd  Chi_s  Vac  he]  for  UAV2 
7.  Cmd  (11:  15)  =  [pNd  pEd  Chi_s  Vac  he]  for  UAV3 

18  7.  Cmd  (16: 20)  =  [pNd  pEd  Chi_s  Vac  he]  for  UAV4 

7»Set  up  target  and  desired  paths 
pN_d  =  200 ; 
pE_d  =  200 ; 

23  pN_i=  [300  250  200  0]  ’ ; 
pE_ i =  [0  0  0  50]  ’  ; 
dpN  =  pN_d -pN_i  ; 
dpE  =  pE_d-pE_i  ; 

Chi_s=atan2(dpE,dpN) ; 

28  he  =50 ; 

7.  Find  commanded  velocities 
kV= . 12 ; 

V_nom  =  11 . 75  ; 

33  dpN_u= [yin (3)  yin(6)  yin(9)  yin  (  12) ]  ; -pN_d ; 
dpE_u= [yin  (4)  yin(7)  yin(10)  yin ( 13) ]  ’ -pE_d ; 
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d_ i = [norm ( [ dpE_u ( 1 )  dpN_u(l)]); 

norm ( [dpE_u (2)  dpN_u(2)]); 

norm ( [ dpE_u (3)  dpN_u(3)]); 

38  norm ( [ dpE_u (4)  dpN_u(4)])]; 

d_mean=mean ( d_ i ) ; 

Vgc = V_nom -kV* ( d_mean - d_i ) ; 

Y=[yin(5)  yin(8)  yin(ll)  yin(14)]’; 

Vw=  [yin ( 1 )  yin (2)  ; yin ( 1)  yin  (2)  ; yin (1)  y in ( 2 )  ; y in ( 1 )  yin(2)]; 

43  Vac_v= [ Vgc . *  cos ( Y)  Vgc . *  sin ( Y) ] - Vw ; 

Vac=[norm(Vac_v(l ,  :))  ;norm(Vac_v(2,  :) )  ; norm (Vac_v(3,  :))  ;norm(Vac_v(4,  :))]  ; 
for  ii=l : 4 

if  Vac ( ii ) >  21.75 
Vac (ii) =21 . 75 ; 

48  end 

end 

7,  Create  output 

Cmd(l:5)  =  [pN_d  pE_d  Chi_s(l)  Vac(l)  he]; 

53  Cmd(6:10)  =  [pN_d  pE_d  Chi_s(2)  Vac (2)  he]; 

Cmd (11  :  15)  =  [pN_d  pE_d  Chi_s(3)  Vac  (3)  he]; 

Cmd(16:20)  =  [pN_d  pE_d  Chi_s  (4)  Vac  (4)  he]; 

Cmd  =  Cmd  ’  ; 


C.3  Mode  3  Controller 

Listing  C.3:  SourceCocle/CC  Controlled. in 

function  Cmd=CC_Controller 1 (yin) 

2  '/.This  function  is  used  in  conjunction  with  Simul ink  model  BATCAM_1  to  model 
'/.the  velocity  feedback  algorithm 
7.  Input  s  : 

7.  yin(l:2)  =  wind  vector  [wNorth  wEast] 

7.  yin(3:5)  =  position  of  UAV1  [pN  pE  Yaw] 

7  7.  yin(6:8)  =  position  of  UAV2  [pN  pE  Yaw] 

1  yin(9:ll)  =  position  of  UAV3  [pN  pE  Yaw] 

7.  yin(12:14)  =  position  of  UAV4  [pN  pE  Yaw] 
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°/o  Out  puts 

7.  [pNd  pEd]  =  desired  location  for  the  UAV 
12  1  Chi_s  =  desired  path  heading 
7.  Vac  =  commanded  airspeed 
7.  he  =  commanded  altitude 

7.  Cmd  (1:5)  =  [pNd  pEd  Chi_s  Vac  he]  for  UAV1 
7.  Cmd  (6:  10)  =  [pNd  pEd  Chi_s  Vac  he]  for  UAV2 
17  7.  Cmd  (11:15)  =  [pNd  pEd  Chi_s  Vac  he]  for  UAV3 

7.  Cmd  (16: 20)  =  [pNd  pEd  Chi_s  Vac  he]  for  UAV4 

7.Set  up  target  and  desired  paths 
pN_d=0 ; 

22  pE_d=0 ; 
he  =50 ; 

7.  Find  commanded  velocities 
V_nom  =  11 . 75  ; 

27  pN_u=[yin(3)  yin(6)  yin(9)  yin(12)]’; 
pE_u=[yin(4)  yin(7)  yin(10)  yin(13)]’; 
pNpE=[pN_u  pE_u] ; 
pNpE_tgt = [pN_d  pE_d]; 

Vgc  =  Maintain0rbit2 (pNpE ,  pNpE_tgt  , V_nom  ,  [ . 5  .1]); 

32  Y=[yin(5)  yin(8)  yin(ll)  yin(14)]’; 

Vw= [yin ( 1 )  yin (2) ; yin ( 1)  yin (2) ; yin (1)  y in ( 2 ) ; y in ( 1 )  yin(2)]; 

Vac_v=[Vgc ’ .  *  cos ( Y)  Vgc ’ ,*sin(Y)]-Vw; 

Vac=[norm(Vac_v(l , :)) ;norm(Vac_v(2, :)) ; norm (Vac_? (3 , :)) ;norm(Vac_v(4, :))] 
for  ii=l : 4 

37  if  Vac ( ii ) >  21.75 

Vac (ii) =21 . 75 ; 

end 

end 

7. Find  Desired  path  heading  Chi_s 
42  i  i  =  1  :  4  ; 

Ang_d=-(ii-l) *pi/2; 

Chi_s=atan2 (pE_u , pN_u) -pi/2; 
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'/.Create  output 

47  Cmd (1:5)  =  [pN_d  pE_d  Chi_s(l)  Vac(l)  he]; 
Cmd (6 : 10)  =  [pN_d  pE_d  Chi_s(2)  Vac (2)  he]; 
Cmd (11 : 15)  =  [pN_d  pE_d  Chi_s  (3)  Vac (3)  he] 
Cmd(16:20)  =  [pN_d  pE_d  Chi_s  (4)  Vac (4)  he] 
Cmd  =  Cmd  ’  ; 
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